Divided Differences: Reply to Hui & Iverson

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

I should like to thank Roger Hui and Ken Iverson for contributing the elegant solutions for divided differences, Choleski and QR decomposition in the January Vector [1]. I would make the point that the MWY solutions reflect my empirical observation that most mathematicians, faced with a programming task, find it more comforting and secure to work at the cell level, rather than to "think array". I don’t applaud this, but it does seem to be a fact that mathematicians, like teenagers, do not always do what they are told is good for them! Cell-level working is also the style of most numerical analysis algorithms given in the older classic text-books.

I used the Choleski algorithm because I wanted a simple, but non-trivial illustration of an algorithm which demanded a loop. Clearly I should have looked further for my example! I can only plead that my upbringing on the books of the pre-computer era has shielded me from the Choleski and QR algorithms in J which Ken and Roger give. These are built on results which are so delightful that I consider it worth while "drawing" them below in ways which make the recursive algorithms seem almost inevitable! C(X) means the Choleski matrix of the square symmetrical matrix X, and is the extension of the idea of "square root" to a matrix. Dashes represent transpose. The matrices are divided into halves either by both rows and columns in the case of Choleski, or by columns in the case of QR. Where the number of columns is odd, the extra one is to the left of the dividing vertical bar.

   X | Y                C(X)         |        0
C(---+---)  =   ---------------------+---------------------
   Y'| Z         Y'X⁻¹ C(X)          |  C(Z - Y'X⁻¹ Y)

                ┌───────┬───────┐   ┌───────┬──────────┐
QR(A0|A1)  =    │ Q(A0) │ Q(G)  │   │ R(A0) │ Q(A0)'A1 │
                │       │       │   ├───────┼──────────┤
                │       │       │   │   0   │   R(G)   │
                └───────┴───────┘   └───────┴──────────┘

where G = A1 - Q(A0).Q(A0)′.A1.

I hope that this further note may encourage others to revisit the J algorithms.

Reference¶

[1] Hui, Roger KW and Iverson, Kenneth, A Note on Programming Style, Vector, Volume 12 no. 3, pp.117-121, January 1996