qr.lisp 11.3 KB
Newer Older
liam's avatar
liam committed
1
2
;; QR decomposition
;; Liam Healy 2008-02-17 11:05:20EST qr.lisp
3
;; Time-stamp: <2009-12-22 22:34:28EST qr.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
liam's avatar
Add QR.    
liam committed
5
6
7

(in-package :gsl)

8
9
;;; /usr/include/gsl/gsl_linalg.h

liam's avatar
liam committed
10
11
12
13
14
;;; FDL
;;; A general rectangular M-by-N matrix A has a
;;; QR decomposition into the product of an orthogonal
;;; M-by-M square matrix Q (where Q^T Q = I) and
;;; an M-by-N right-triangular matrix R, A = Q R.
liam's avatar
Add QR.    
liam committed
15

liam's avatar
liam committed
16
17
18
19
20
21
;;; This decomposition can be used to convert the linear system A x = b
;;; into the triangular system R x = Q^T b, which can be solved by
;;; back-substitution. Another use of the QR decomposition is to
;;; compute an orthonormal basis for a set of vectors. The first N
;;; columns of Q form an orthonormal basis for the range of A,
;;; ran(A), when A has full column rank.
liam's avatar
Add QR.    
liam committed
22

23
24
25
26
(defmfun QR-decomposition
    (A
     &optional
     (tau (make-marray 'double-float :dimensions (min (dim0 A) (dim1 A)))))
liam's avatar
Add QR.    
liam committed
27
  "gsl_linalg_QR_decomp"
Liam Healy's avatar
Liam Healy committed
28
29
30
  (((mpointer A) :pointer) ((mpointer tau) :pointer))
  :inputs (A)
  :outputs (A tau)
liam's avatar
liam committed
31
32
33
  :documentation 			; FDL
  "Factorize the M-by-N matrix A into the QR decomposition A = Q R.
   On output the diagonal and
liam's avatar
Add QR.    
liam committed
34
   upper triangular part of the input matrix contain the matrix
liam's avatar
liam committed
35
36
37
38
39
   R.  The vector tau and the columns of the lower triangular
   part of the matrix A contain the Householder coefficients and
   Householder vectors which encode the orthogonal matrix Q.  The
   vector tau must be of length k=min(M,N). The matrix
   Q is related to these components by, Q = Q_k ... Q_2 Q_1
Liam Healy's avatar
Liam Healy committed
40
   where Q_i = I - tau_i v_i v_i^T and v_i is the
liam's avatar
liam committed
41
42
   Householder vector v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)).
   This is the same storage scheme as used by lapack.
liam's avatar
Add QR.    
liam committed
43
44

   The algorithm used to perform the decomposition is Householder QR (Golub
liam's avatar
liam committed
45
   & Van Loan, Matrix Computations, Algorithm 5.2.1).")
liam's avatar
Add QR.    
liam committed
46

47
48
49
50
51
52
53
54
55
56
57
58
(defmfun QR-solve
    (QR tau b &optional x-spec
       &aux
       (x (if (eq x-spec t)
	      (make-marray 'double-float :dimensions (dimensions b))
	      x-spec)))
  ("gsl_linalg_QR_svx" "gsl_linalg_QR_solve")
  ((((mpointer QR) :pointer) ((mpointer tau) :pointer)
    ((mpointer b) :pointer))
   (((mpointer QR) :pointer) ((mpointer tau) :pointer)
    ((mpointer b) :pointer) ((mpointer x) :pointer)))
  :inputs (QR tau b x)
Liam Healy's avatar
Liam Healy committed
59
  :outputs (x)
60
  :return ((or x b))
liam's avatar
liam committed
61
  :documentation			; FDL
62
63
64
65
66
67
68
69
70
71
  "Solve the square system A x = b using the QR decomposition of A
   into (QR, tau) given by QR-decomp. The least-squares solution for
   rectangular systems can be found using QR-lssolve.  If x-spec is
   NIL (default), the solution will replace b.  If x-spec is T, then
   an array will be created and the solution returned in it.  If
   x-spec is a marray, the solution will be returned in it.  If x-spec
   is non-NIL, on output the solution is stored in x and b is not
   modified.  The solution is returned from the function call.")

(defmfun QR-solve-least-squares
72
73
74
75
    (QR tau b
	&optional
	 (x (make-marray 'double-float :dimensions (dim1 QR)))
	 (residual (make-marray 'double-float :dimensions (dim0 QR))))
liam's avatar
Add QR.    
liam committed
76
  "gsl_linalg_QR_lssolve"
Liam Healy's avatar
Liam Healy committed
77
78
79
80
81
  (((mpointer QR) :pointer) ((mpointer tau) :pointer)
   ((mpointer b) :pointer) ((mpointer x) :pointer)
   ((mpointer residual) :pointer))
  :inputs (QR tau b)
  :outputs (x residual)
liam's avatar
liam committed
82
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
83
84
85
86
87
88
89
90
  "The least squares solution to the overdetermined system A x = b
   where the matrix A has more rows than columns.  The least squares
   solution minimizes the Euclidean norm of the residual, ||Ax -
   b||.The routine uses the QR decomposition of A into (QR, tau) given
   by #'QR-decomposition.  The solution is returned in x.  The
   residual is computed as a by-product and stored in residual.")

(defmfun QR-QTvector (QR tau v)
91
  "gsl_linalg_QR_QTvec"
Liam Healy's avatar
Liam Healy committed
92
93
94
95
  (((mpointer QR) :pointer) ((mpointer tau) :pointer)
   ((mpointer v) :pointer))
  :inputs (QR tau)
  :outputs (v)
liam's avatar
liam committed
96
97
98
99
100
101
102
  :documentation			; FDL
  "Apply the matrix Q^T encoded in the decomposition
   (QR, tau) to the vector v, storing the result Q^T v in v.
   The matrix multiplication is carried out directly using
   the encoding of the Householder vectors without needing to form the full
   matrix Q^T.")

Liam Healy's avatar
Liam Healy committed
103
(defmfun QR-Qvector (QR tau v)
104
  "gsl_linalg_QR_Qvec"
Liam Healy's avatar
Liam Healy committed
105
106
107
108
  (((mpointer QR) :pointer) ((mpointer tau) :pointer)
   ((mpointer v) :pointer))
  :inputs (QR tau)
  :outputs (v)
liam's avatar
liam committed
109
110
111
112
  :documentation			; FDL
  "Apply the matrix Q encoded in the decomposition
   (QR, tau) to the vector v, storing the result Q v in v.
   The matrix multiplication is carried out directly using
liam's avatar
Add QR.    
liam committed
113
   the encoding of the Householder vectors without needing to form the full
liam's avatar
liam committed
114
   matrix Q.")
liam's avatar
Add QR.    
liam committed
115

116
117
118
119
120
121
122
123
124
125
126
127
(defmfun QR-Rsolve
    (QR b &optional x-spec
       &aux
       (x (if (eq x-spec t)
	      (make-marray 'double-float :dimensions (dimensions b))
	      x-spec)))
  ("gsl_linalg_QR_Rsvx" "gsl_linalg_QR_Rsolve")
  ((((mpointer QR) :pointer) ((mpointer b) :pointer))
   (((mpointer QR) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer)))
  :inputs (QR b x)
  :outputs (b x)
  :return ((or x b))
liam's avatar
liam committed
128
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
129
  "Solve the triangular system R x = b for x.  It may be useful if the
130
131
132
133
134
135
   product b' = Q^T b has already been computed using QR-QTvec.  If
   x-spec is NIL (default), the solution will replace b.  If x-spec is
   T, then an array will be created and the solution returned in it.
   If x-spec is a marray, the solution will be returned in it.  If
   x-spec is non-NIL, on output the solution is stored in x and b is
   not modified.  The solution is returned from the function call.")
liam's avatar
Add QR.    
liam committed
136

137
138
139
140
141
(defmfun QR-unpack
    (QR tau
     &optional
     (Q (make-marray 'double-float :dimensions (list (dim0 QR) (dim0 QR))))
     (R (make-marray 'double-float :dimensions (dimensions QR))))
liam's avatar
Add QR.    
liam committed
142
  "gsl_linalg_QR_unpack"
Liam Healy's avatar
Liam Healy committed
143
144
145
146
  (((mpointer QR) :pointer) ((mpointer tau) :pointer)
   ((mpointer Q) :pointer) ((mpointer R) :pointer))
  :inputs (QR tau)
  :outputs (Q R)
liam's avatar
liam committed
147
148
149
150
  :documentation			; FDL
  "Unpack the encoded QR decomposition
  (QR, tau) into the matrices Q and R where
  Q is M-by-M and R is M-by-N.")
liam's avatar
Add QR.    
liam committed
151

152
153
(defmfun QR-QRsolve
    (Q R b &optional (x (make-marray 'double-float :dimensions (dim0 b))))
liam's avatar
Add QR.    
liam committed
154
  "gsl_linalg_QR_QRsolve"
Liam Healy's avatar
Liam Healy committed
155
156
157
158
  (((mpointer Q) :pointer) ((mpointer R) :pointer)
   ((mpointer b) :pointer) ((mpointer x) :pointer))
  :inputs (Q R b)
  :outputs (x)
liam's avatar
liam committed
159
160
161
162
  :documentation			; FDL
  "Solves the system R x = Q^T b for x.  It can
  be used when the QR decomposition of a matrix is available in
  unpacked form as Q, R).")
liam's avatar
Add QR.    
liam committed
163

liam's avatar
liam committed
164
(defmfun QR-update (Q R w v)
liam's avatar
Add QR.    
liam committed
165
  "gsl_linalg_QR_update"
Liam Healy's avatar
Liam Healy committed
166
167
168
169
  (((mpointer Q) :pointer) ((mpointer R) :pointer)
   ((mpointer w) :pointer) ((mpointer v) :pointer))
  :inputs (Q R w v)
  :outputs (w Q R)
liam's avatar
liam committed
170
171
172
173
174
175
176
177
  :return (Q R)
  :documentation			; FDL
  "Perform a rank-1 update w v^T of the QR
  decomposition (Q, R). The update is given by Q'R' = Q R + w v^T
  where the output matrices Q' and R' are also
  orthogonal and right triangular. Note that w is destroyed by the
  update.")

178
179
180
181
182
183
184
185
186
187
(defmfun R-solve
    (R b &optional x-spec
       &aux
       (x (if (eq x-spec t)
	      (make-marray 'double-float :dimensions (dimensions b))
	      x-spec)))
  ("gsl_linalg_R_svx" "gsl_linalg_R_solve")
  ((((mpointer R) :pointer) ((mpointer b) :pointer))
   (((mpointer R) :pointer) ((mpointer b) :pointer) ((mpointer x) :pointer)))
  :inputs (R b x)
Liam Healy's avatar
Liam Healy committed
188
  :outputs (x)
189
  :return ((or x b))
liam's avatar
liam committed
190
  :documentation			; FDL
191
192
193
194
195
196
197
198
  "Solve the triangular system R x = b in-place. On input x should
  contain the right-hand side b, which is replaced by the solution on
  output.  If x-spec is NIL (default), the solution will replace b.
  If x-spec is T, then an array will be created and the solution
  returned in it.  If x-spec is a marray, the solution will be
  returned in it.  If x-spec is non-NIL, on output the solution is
  stored in x and b is not modified.  The solution is returned from
  the function call.")
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245

;;; Examples and unit test, from linalg/test.c

(defun test-qr-solve-dim (matrix)
  "Solve the linear equation using QR with the supplied matrix and
   a right-hand side vector which is the reciprocal of one more than
   the index."
  (let ((dim (dim0 matrix)))
    (multiple-value-bind (QR tau)
	(QR-decomposition (copy matrix))
      (QR-solve QR tau (create-rhs-vector dim) T))))

(defun test-qr-qrsolve-dim (matrix)
  "Solve the linear equation using QR with the supplied matrix and
   a right-hand side vector which is the reciprocal of one more than
   the index."
  (let ((dim (dim0 matrix)))
    (multiple-value-bind (QR tau)
	(QR-decomposition (copy matrix))
      (multiple-value-bind (Q R)
	  (QR-unpack QR tau)
	(QR-QRsolve Q R (create-rhs-vector dim))))))

(defun test-qr-lssolve-dim (matrix)
  "Solve the linear equation using QR least squares with the supplied
   matrix and a right-hand side vector which is the reciprocal of one
   more than the index.  Returns the solution and the residual."
  (let ((dim (dim0 matrix)))
    (multiple-value-bind (QR tau)
	(QR-decomposition (copy matrix))
      ;; Residual not checked.
      (QR-solve-least-squares QR tau (create-rhs-vector dim)))))

(defun test-qr-decomp-dim (matrix)
  "Solve the QR decomposition with the supplied
   matrix and a right-hand side vector which is the reciprocal of one
   more than the index."
  (multiple-value-bind (QR tau)
      (QR-decomposition (copy matrix))
    (multiple-value-bind (Q R)
	(QR-unpack QR tau)
      (matrix-product Q R))))

(defun test-qr-update-dim (matrix)
  "Test QR rank-1 update; this should return a matrix with all
   elements near zero."
  (let* ((dim0 (dim0 matrix)) (dim1 (dim1 matrix))
246
247
	 (u (create-matrix (lambda (i) (sin (1+ i))) dim0 nil))
	 (v (create-matrix
248
	     (lambda (i) (+ (cos (+ 2 i)) (sin (+ 3 (expt i 2)))))
249
	     dim1 nil))
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
	 (qr1
	  (create-matrix
	   (lambda (i j) (+ (maref matrix i j) (* (maref u i) (maref v j))))
	   dim0 dim1))
	 (qr2 (copy matrix))
	 (w (make-marray 'double-float :dimensions dim0)))
    (multiple-value-bind (QR2 tau)
	(QR-decomposition qr2)
      (multiple-value-bind (Q2 R2)
	  (QR-unpack QR2 tau)
	;; compute w = Q^T u
	(matrix-product Q2 u w 1.0d0 0.0d0 :trans)
	(QR-update Q2 R2 w v)
	(matrix-product Q2 R2 qr2 1.0d0 0.0d0)
	(elt- qr1 qr2)))))

(save-test qr
 ;; test_QR_solve
 (test-qr-solve-dim *hilb2*)
 (test-qr-solve-dim *hilb3*)
 (test-qr-solve-dim *hilb4*)
 (test-qr-solve-dim *hilb12*)
 (test-qr-solve-dim *vander2*)
 (test-qr-solve-dim *vander3*)
 (test-qr-solve-dim *vander4*)
 (test-qr-solve-dim *vander12*)
 ;; test_QR_QRsolve
 (test-qr-qrsolve-dim *hilb2*)
 (test-qr-qrsolve-dim *hilb3*)
 (test-qr-qrsolve-dim *hilb4*)
 (test-qr-qrsolve-dim *hilb12*)
 (test-qr-qrsolve-dim *vander2*)
 (test-qr-qrsolve-dim *vander3*)
 (test-qr-qrsolve-dim *vander4*)
 (test-qr-qrsolve-dim *vander12*)
 ;; test_QR_lssolve
 (test-qr-lssolve-dim *m53*)
 (test-qr-lssolve-dim *hilb2*)
 (test-qr-lssolve-dim *hilb3*)
 (test-qr-lssolve-dim *hilb4*)
 (test-qr-lssolve-dim *hilb12*)
 (test-qr-lssolve-dim *vander2*)
 (test-qr-lssolve-dim *vander3*)
 (test-qr-lssolve-dim *vander4*)
 (test-qr-lssolve-dim *vander12*)
 ;; test_QR_decomp
 (test-qr-decomp-dim *m35*)
 (test-qr-decomp-dim *m53*)
 (test-qr-decomp-dim *hilb2*)
 (test-qr-decomp-dim *hilb3*)
 (test-qr-decomp-dim *hilb4*)
 (test-qr-decomp-dim *hilb12*)
 (test-qr-decomp-dim *vander2*)
 (test-qr-decomp-dim *vander3*)
 (test-qr-decomp-dim *vander4*)
 (test-qr-decomp-dim *vander12*)
 ;; test_QR_update
 (test-qr-update-dim *m35*)
 (test-qr-update-dim *m53*)
 (test-qr-update-dim *hilb2*)
 (test-qr-update-dim *hilb3*)
 (test-qr-update-dim *hilb4*)
 (test-qr-update-dim *hilb12*)
 (test-qr-update-dim *vander2*)
 (test-qr-update-dim *vander3*)
 (test-qr-update-dim *vander4*)
 (test-qr-update-dim *vander12*))