A Note on Programming Style 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 WMY 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 Jottings 7, the WMY program can be made even more Pascallike 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 2by2 matrix of matrices and substituting appropriate matrix operations for scalar ones in the Choleski method for a 2by2 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 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 2column matrix of matrices and substituting matrix operations for scalar and vector ones in the GramSchmidt 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 2element 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 1References
Originally appeared in Vector, Volume 12, Number 3, 199601.
