Matrix Decomposition

by Norman Thomson

mathematics and statistics programming techniques

Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 31

J contains, as a primitive, 128!:0, the so-called QR decomposition of a square matrix. This note describes how to obtain this in APL, and also how to obtain the LU decomposition where L and U stand for Lower and Upper triangular matrices, in the latter case with 1s on the diagonal. L and U can if required be merged compactly into a matrix of the same shape as the orginal matrix by overlaying all those cells which, by definition, must be either 0 or 1.

QR decomposition is also known as Gram-Schmidt orthogonalisation. The comments on the functions should make the procedures clear, particularly if used alongside the sort of algebraic descriptions of the algorithms which can be found in many texts on Numerical Analysis such as "Matrix Computations" by Golub and Van Loan, or the "Numerical Recipes" series. Perhaps surprisingly, it is in the nature of the algorithms that they are not particularly helped by the availability of the array operations of APL.

    ∇ Z←QR R;T;U;I;N
[1]  ⍝ R is square.  R≡⊃+.×/Z, cols of ↑Z are orthonormal, 2⊃Z is upper ∆r
[2]   N←↑⍴R ⋄ Z←2⍴⊂(N,0)⍴I←0             ⍝ initialize, Z is built col by col
[3]  L1:→(N<I←I+1)/0 ⋄ T←R[;I]           ⍝ begin loop, NORM is sqrt +/x*2
[4]   Z←1E¯6 CLEAN Z,¨(⊂U÷NORM U),⊂N↑T,NORM U←+/(1,-T←T+.×↑Z)×[2]T,↑Z ⋄ →L1
    ∇

    ∇ Z←NORM R
[1]   Z←(+/R*2)*0.5                      ⍝ Euclidean norm
    ∇

    ∇ Z←L CLEAN R
[1]   Z←R×L≤|R                           ⍝ sets small numbers to zero
    ∇

    ∇ Z←LU R;I;N;T;K;U
[1]  ⍝ ↑Z is lower ∆r, 2⊃Z is upper ∆r with diag 1s.  R≡⊃+.×/LU Z
[2]   Z←R ⋄ U←⍳N←↑⍴R ⋄ I←K←0                   ⍝ initialize for loop[3-6]
[3]  L1:→(N≤I←I+1)/L3 ⋄ T←I↓U ⋄ →(I=1)/L2     ⍝ first row is special case
[4]   K←Z[,I;⍳I-1]+.×Z[⍳I-1;T]                 ⍝ K is for use in [5]
[5]  L2:Z[I;T]←(R[,I;T]-K)÷Z[I;I]              ⍝ set Ith row to rt of diag
[6]   Z[T;I+1]←R[T;I+1]-Z[T;⍳I]+.×Z[⍳I;I+1] ⋄ →L1  ⍝ diag+below of (I+1)th col
[7]  L3:Z←((U∘.≥U)×Z)((U∘.=U)+(U∘.<U)×Z)       ⍝ omit [7] for compact form
    ∇

      M
 8 2 1 3
 1 5 1 8
10 3 2 4
 9 7 2 7

      QR M
0.5101  ¯0.2921 ¯0.5667  0.5774   15.68 7.268 2.997  8.607
0.06376  0.7761  0.2456  0.5774    0    5.846 0.8929 7.602
0.6376  ¯0.2796  0.7179  0         0    0     0.4723 0.8879
0.5738   0.484  ¯0.3212 ¯0.5774    0    0     0      2.309

      LU M
 8 0    0       0    1 0.25 0.125   0.375
 1 4.75 0       0    0 1    0.1842  1.605
10 0.5  0.6579  0    0 0    1      ¯0.84
 9 4.75 0      ¯4    0 0    0       1