QR decomposition

Any rectangular $m\times n$ matrix $A$ can be decomposed to a product of a orthogonal matrix $Q$ and a upper (right) triangular matrix $R$, as described in QR decomposition.

There are several methods to compute $Q$ and $R$, but in this task, the method of Householder reflections should be used. Multiplying a given vector $a$, for example the first column of matrix $A$, with the Householder matrix $H$, which is given as

$H=I-{\frac {2}{u^{T}u}}uu^{T}$

reflects $a$ about a plane given by its normal vector $u$. When the normal vector of the plane $u$ is given as

$u=a-\|a\|_{2}\;e_{1}$

then the transformation reflects $a$ onto the first standard basis vector

$e_{1}=[1\;0\;0\;...]^{T}$

which means that all entries but the first become zero. To avoid numerical cancellation errors, we should take the opposite sign of $a_{1}$:

$u=a+{\textrm {sign}}(a_{1})\|a\|_{2}\;e_{1}$

and normalize with respect to the first element:

$v={\frac {u}{u_{1}}}$

Applying $H$ on $a$ then gives

$H\;a=-{\textrm {sign}}(a_{1})\;\|a\|_{2}\;e_{1}$

and applying $H$ on the matrix $A$ zeroes all subdiagonal elements of the first column:

$H_{1}\;A={\begin{pmatrix}r_{11}&r_{12}&r_{13}\\0&*&*\\0&*&*\end{pmatrix}}$

In the second step, the second column of $A$, we want to zero all elements but the first two, which means that we have to calculate $H$ with the first column of the submatrix (denoted *), not on the whole second column of $A$.

To get $H_{2}$, we then embed the new $H$ into an $m\times n$ identity:

$H_{2}={\begin{pmatrix}1&0&0\\0&H&\\0&&\end{pmatrix}}$

This is how we can, column by column, remove all subdiagonal elements of $A$ and thus transform it into $R$.

$H_{n}\;...\;H_{3}H_{2}H_{1}A=R$

The product of all the Householder matrices $H$, for every column, in reverse order, will then yield the orthogonal matrix $Q$.

$H_{1}H_{2}H_{3}\;...\;H_{n}=Q$

The QR decomposition should then be used to solve linear least squares (Multiple regression) problems $Ax=b$ by solving

$R\;x=Q^{T}\;b$

When $R$ is not square, i.e. $m>n$ we have to cut off the $m-n$ zero padded bottom rows.

$R={\begin{pmatrix}R_{1}\\0\end{pmatrix}}$

and the same for the RHS:

$Q^{T}\;b={\begin{pmatrix}q_{1}\\q_{2}\end{pmatrix}}$

Finally, solve the square upper triangular system by back substitution:

$R_{1}\;x=q_{1}$

Demonstrate the QR decomposition on the example matrix from the Wikipedia article:

$A={\begin{pmatrix}12&-51&4\\6&167&-68\\-4&24&-41\end{pmatrix}}$

and the usage for linear least squares problems on the example from Polynomial_regression.

Common Lisp

Uese the routines m+, m-, .* from Element-wise_operations, mmul from Matrix multiplication, mtp from Matrix transposition.

Helper functions: <lang lisp>(defun sign (x)

 (if (zerop x)
     x
     (/ x (abs x))))

(defun norm (x)

 (let ((len (car (array-dimensions x))))
   (sqrt (loop for i from 0 to (1- len) sum (expt (aref x i 0) 2)))))

(defun make-unit-vector (dim)

 (let ((vec (make-array `(,dim ,1) :initial-element 0.0d0)))
   (setf (aref vec 0 0) 1.0d0)
   vec))

(defun array-range (A ma mb na nb)

 (let* ((m (car (array-dimensions A)))
        (n (car (array-dimensions A)))
        (mm (1+ (- mb ma)))
        (nn (1+ (- nb na)))
        (B (make-array `(,mm ,nn) :initial-element 0.0d0)))

   (loop for i from 0 to (1- mm) do
        (loop for j from 0 to (1- nn) do
             (setf (aref B i j)
                   (aref A (+ ma i) (+ na j)))))
   B))

(defun rows (A) (car (array-dimensions A))) (defun cols (A) (cadr (array-dimensions A))) (defun mcol (A n) (array-range A 0 (1- (rows A)) n n)) (defun mrow (A n) (array-range A n n 0 (1- (cols A))))

(defun array-embed (A B row col)

 (let* ((ma (rows A))
        (na (cols A))
        (mb (rows B))
        (nb (cols B))
        (C  (make-array `(,ma ,na) :initial-element 0.0d0)))

   (loop for i from 0 to (1- ma) do
        (loop for j from 0 to (1- na) do
             (setf (aref C i j) (aref A i j))))

   (loop for i from 0 to (1- mb) do
        (loop for j from 0 to (1- nb) do
             (setf (aref C (+ row i) (+ col j))
                   (aref B i j))))

   C))

</lang>

Main routines: <lang lisp> (defun make-householder (a)

 (let* ((m    (car (array-dimensions a)))
        (s    (sign (aref a 0 0)))
        (e    (make-unit-vector m))
        (u    (m+ a (.* (* (norm a) s) e)))
        (v    (./ u (aref u 0 0)))
        (beta (/ 2 (mmul (mtp v) v))))
   
   (m- (eye m)
       (.* beta (mmul v (mtp v))))))

(defun qr (A)

 (let* ((m (car  (array-dimensions A)))
        (n (cadr (array-dimensions A)))
        (Q (eye m)))

   ;; Work on n columns of A.
   (loop for i from 0 to (if (= m n) (- n 2) (- n 1)) do

        ;; Select the i-th submatrix. For i=0 this means the original matrix A.
        (let* ((B (array-range A i (1- m) i (1- n)))
               ;; Take the first column of the current submatrix B.
               (x (mcol B 0))
               ;; Create the Householder matrix for the column and embed it into an mxm identity.
               (H (array-embed (eye m) (make-householder x) i i)))

          ;; The product of all H matrices from the right hand side is the orthogonal matrix Q.
          (setf Q (mmul Q H))

          ;; The product of all H matrices with A from the LHS is the upper triangular matrix R.
          (setf A (mmul H A))))

   ;; Return Q and R.
   (values Q A)))

</lang>

Example 1:

<lang lisp>(qr #2A((12 -51 4) (6 167 -68) (-4 24 -41)))

1. 2A((-0.85 0.39 0.33)

   (-0.42 -0.90 -0.03)
   ( 0.28 -0.17  0.94))

1. 2A((-14.0 -21.0 14.0)

   (  0.0 -175.0  70.0)
   (  0.0    0.0 -35.0))</lang>

Example 2, Polynomial regression:

<lang lisp>(defun polyfit (x y n)

 (let* ((m (cadr (array-dimensions x)))
        (A (make-array `(,m ,(+ n 1)) :initial-element 0)))
   (loop for i from 0 to (- m 1) do
         (loop for j from 0 to n do
               (setf (aref A i j)
                     (expt (aref x 0 i) j))))
   (lsqr A (mtp y))))

Solve a linear least squares problem by QR decomposition.

(defun lsqr (A b)

 (multiple-value-bind (Q R) (qr A)
   (let* ((n (cadr (array-dimensions R))))
     (solve-upper-triangular (array-range R                0 (- n 1) 0 (- n 1))
                             (array-range (mmul (mtp Q) b) 0 (- n 1) 0 0)))))

Solve an upper triangular system by back substitution.

(defun solve-upper-triangular (R b)

 (let* ((n (cadr (array-dimensions R)))
        (x (make-array `(,n 1) :initial-element 0.0d0)))

   (loop for k from (- n 1) downto 0
      do (setf (aref x k 0)
               (/ (- (aref b k 0)
                     (loop for j from (+ k 1) to (- n 1)
                        sum (* (aref R k j)
                               (aref x j 0))))
                  (aref R k k))))
   x))</lang>

<lang lisp>;; Finally use the data: (let ((x #2A((0 1 2 3 4 5 6 7 8 9 10)))

     (y #2A((1 6 17 34 57 86 121 162 209 262 321))))
   (polyfit x y 2))

1. 2A((0.999999966345088) (2.000000015144699) (2.99999999879804))</lang>