A Note on Programming Style

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

We were pleased to see the uses of J in the October Education Vector, particularly the Muller, van Woudenberg, and Young program for divided differences and the fact that that it resulted from a mathematics program at Strathclyde. We would like to use it together with contributions from Norman Thomson in the same issue to discuss certain matters of programming style: although MWY followed "the style of a Pascal program, which is a natural way …", the use of vectors and matrices would be more natural to mathematicians, and would better lend itself to analysis.

Divided Differences¶

As shown by the use of control structures in Thomson’s J-ottings 7, the MWY program can be made even more Pascal-like as follows:

dd0 =: 4 : 0
 i=.0
 n=.#x.
 z=.(1<.n){.d=.y.
 while. n>i=.1+i do.
  dx=.(-i)}.(i|.x.)-x.
  dy=.2 --/\ d
  d =.dy%dx
  z =.z,{.d
 end.
)

   x=.0 1 5 8
   y=.4 6 18 6
   x dd0 y
4 2 0.2 _0.15

The expression for dx computes differences in x that are i apart, by taking the difference between rotating x by i and x itself and discarding the last i entries. The expression for dy divides d into overlapping windows of size 2 and applies --/ (insert subtract from) to each window.

The computation can also be expressed tacitly, an expression that facilitates investigation of the internal workings of the algorithm:

dix=: -@[ }. |. - ]
dx =: >:@-&# dix [
dy =: 2: --/\ ]
dd =: (dy % dx)^:(i.@#@])

   x dd y              {."1 x dd y
    4  6 18 6      4 2 0.2 _0.15
    2  3 _4 0
  0.2 _1  0 0
_0.15  0  0 0

   x dx y              1 dix x
1 4 3                    1 4 3

   x dy y              2 dix x
2 12 _12           5 7

                       3 dix x
                   8

   f=.dy % dx

   x f y
2 3 _4

   x f x f y
0.2 _1

   x f x f x f y
_0.15

   x f^:3 y
_0.15

   x f^:0 1 2 3 y
    4  6 18 6
    2  3 _4 0
  0.2 _1  0 0
_0.15  0  0 0

Choleski Decomposition¶

MWY and Thomson presented programs for the Choleski decomposition. We present here a recursive version. Given a positive definite matrix A, the Choleski method computes a lower triangular matrix L such that A-:L x h L, where x=:+/ .* is the matrix product and h=:+@|: is the conjugate transpose. A recursive solution reveals itself by considering A as a 2-by-2 matrix of matrices and substituting appropriate matrix operations for scalar ones in the Choleski method for a 2-by-2 matrix of scalars.

x=: +/ . *                        matrix product
h=: +@|:                          conjugate transpose

Choleski =: 3 : 0
 n=.#A=.y.
 if. 1>:n do.
  13!:8(,(A=|A)>0=A)}.12         check for positive definite
  %:A
 else.
  p=.>.n%2 [ q=.<.n%2
  X=.(p,p){.A [ Y=.(p,-q){.A [ Z=.(-q,q){.A
  L0=.Choleski X
  L1=.Choleski Z-(T=.(h Y) x %.X) x Y
  L0,(T x L0),.L1
 end.
)

   A                              A random positive definite matrix
  33    7j_8 7j_10     3j_4
 7j8      28   2j4 _10j_11
7j10    2j_4    22     3j3
 3j4 _10j11   3j_3      16

   L=. Choleski A

   $L
4 4

   A -: L x h L                   A true decomposition
1

   *./,(>:/~i.#L) >: 0~:L         L is lower triangular
1

   [ B=. 3 3$ 1 4 6  4 0 3  6 3 2
1 4 6
4 0 3
6 3 2

   Choleski B                     B is not positive definite
|assertion failure
|       13!:8(,(A=|A)>0=A)}.12

QR Decomposition¶

Thomson’s Technical Note on Matrix Decomposition presented the QR and LU decompositions. Both are amenable to the recursive approach used above for the Choleski decomposition, with array operations playing a key role. Here, we present a solution for the QR decomposition; a recursive array LU decomposition can be found in Aho, Hopcroft, and Ullman [1974, Section 6.4].

Given a matrix A where >:/$A, the QR decomposition produces an orthogonal matrix Q and a square upper triangular matrix R such that A-:Q x 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.

x=: +/ . *                        matrix product
h=: +@|:                          conjugate transpose

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

   [ A=. j./_8+?2 7 4$20          A random matrix
 _6j6   7j10   1j7 2j_3
_4j_8   _8j6  5j_2  5j4
 10j7  _1j11  2j_1 8j_4
_8j11   _7j6   2j7  5j5
_8j_7   _1j4  _7j9 0j_3
    5    3j7  10j1 8j_4
 2j_3  _7j_1  5j_5  0j1

   $QR A                          Result is a 2-element boxed vector
2

   'QR'=.QR A                     Assign first box to Q, second to R

   $Q
7 4

   $R
4 4

   A -: Q x R                     A true decomposition
1

   >./,|(=i.{:$Q)-(h Q) x Q       Q is orthogonal
8.51083e_16

   *./,(<:/~i.#R) >: 0~:R         R is upper triangular
1

References¶

[1] Aho, Alfred V., John Hopcroft, and Jeffrey D. Ullman, The Design and Analysis of Computer Algorithms, Addison-Wesley, 1974.

[2] Muller, Antje, Tineke van Woudenberg, and Alistair Young, Two Numerical Algorithms in J, The Education Vector, Volume 12, Number 2, 1995, pp 26-30.

[3] Strang, Gilbert, Linear Algebra and its Applications, Academic Press, 1976.

[4] Thomson, Norman, J-ottings 7, The Education Vector, Volume 12, Number 2, 1995, pp 21-25

[5] Thomson, Norman, Technical Note on Matrix Decomposition, The Education Vector, Volume 12, Number 2, 1995, pp 31-32.