Essays/QR Decomposition

Given a matrix A where >:/$A , the QR decomposition produces an orthogonal matrix Q and a square upper triangular matrix R such that

mp=: +/ . *  NB. matrix product
h =: +@|:    NB. conjugate transpose
 
A -: Q mp R

A recursive solution reveals itself by considering A as a 2-column matrix of matrices and substituting matrix operations for scalar and vector ones in the Gram-Schmidt method for two vectors.

QR=: 3 : 0
 n=.{:$A=.y
 if. 1>:n do.
  A ((% {.@,) ; ]) %:(h A) mp A
 else.
  m =.>.n%2
  A0=.m{."1 A
  A1=.m}."1 A
  'Q0 R0'=.QR A0
  'Q1 R1'=.QR A1 - Q0 mp T=.(h Q0) mp A1
  (Q0,.Q1);(R0,.T),(-n){."1 R1
 end.
)

The verb is now applied to an example:

   ] A=: j./ _8 + 2 7 4 ?.@$ 20   NB. a random matrix
 _2j8  7j_2 11j_3   4j9
    6    11 _8j_6   9j3
_8j11   6j9 _2j11  10j1
 5j_7 10j_4   3j5   4j1
 10j9  _8j4  2j_4  _6j5
  3j9   8j7  _8j5    _4
 _4j5   3j8  5j_3 _5j_1

   $ QR A                         NB. result is a 2-element boxed vector
2
   'Q R'=: QR A                   NB. assign first box to Q, second to R
   $Q
7 4
   $R
4 4

   A -: Q mp R                    NB. a true decomposition
1
   >./, | (=i.{:$Q) - (h Q) mp Q  NB. Q is orthogonal
2.22045e_16
   (0=R) >: >/~i.#R               NB. R is upper triangular
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1

QRcheck has the same argument and result as QR , but incorporates checks:

QRcheck=: 3 : 0
 A=. y
 assert. >:/$A
 'Q R'=.z=. QR A
 'm n'=. $A
 assert. A -: Q mp R
 assert. (($Q) -: m,n) *. 5e_14 > >./, | (=i.n) - (h Q) mp Q
 assert. (($R) -: n,n) *. (0=R) >: >/~i.n
 z
)

See also

Contributed by RogerHui. Substantially the same text previously appeared as part of A Note on Programming Style in Vector, Volume 12, Number 3, January 1996.