Linear Recurrences and Matrix Powers

From: Roger Hui, 29th February 1996

mathematics and statistics programming techniques

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

I read with enjoyment and admiration accounts of Eugene McDonnell’s excellent adventure in "Heron’s Rule & Integer-Area Triangles" in the January Vector (12.3 pp 133—142). The article discusses Sloane’s sequence 700, s=: 1 2 7 26 97 362 ... with recurrence relation

   A(n) = 4A(n-1) - A(n-2)

where (2*s)+/_1 0 1 are the sides of integer-area triangles. Computing the sequence using the double recursion is found to be extremely slow, and is made faster using generating functions, Taylor’s series, and partial fractions.

There is another fast solution, based on array operations. The recurrence can be written as the matrix equation:

┌      ┐     ┌       ┐ ┌      ┐
│A(n-1)│     │ 0   1 │ │A(n-2)│
│      │  =  │       │ │      │
│A(n)  │     │ _1  4 │ │A(n-1)│
└      ┘     └       ┘ └      ┘

This form makes clear why it’s called a linear recurrence. In J,

   x=: +/ .*                         Matrix multiplication
   M=: 0 1,:_1 4
   A0=: 1 2

   M x A0
2 7
   M x M x A0                           (M x M) x A0
7 26                                 7 26
   M x M x M x A0                       (M x M x M) x A0
26 97                                26 97

   M&x^:(i.5) A0                        (M&x^:(i.5) =i.2) x A0
 1   2                                1   2
 2   7                                2   7
 7  26                                7  26
26  97                               26  97
97 362                               97 362

Thus the n-th element of the sequence obtains by n repeated multiplications by M or, equivalently, by multiplying by the n-th power of M. The n-th power of a matrix can be computed by repeated squaring, with a consequent reduction from O(n) operations to O(log n) operations.

We illustrate the process for M to the 30-th power:

   [ b=: #: 30                       Binary representation of 30
1 1 1 1 0

   b # i.-#b                         Powers corresponding to 1s in b
4 3 2 1

   M x M                             Squaring a matrix
_1  4
_4 15

   x~ M
_1  4
_4 15

   x~^:4 3 2 1 M                     Powers with squaring as the function
_109552575  408855776
_408855776 1525870529

     _2911      10864
    _10864      40545

       _15         56
       _56        209

        _1          4
        _4         15

   (M&x^:16 8 4 2 =i.2) -: x~^:4 3 2 1 M
1

   x/ x~^:4 3 2 1 M                  Product of selected powers
_1.11401e16 4.15753e16
_4.15753e16 1.55161e17

   pow=: 4 : 'x/ x~^:(b#i.-#b=.#:y.) x.'

   M pow 30
_1.11401e16 4.15753e16
_4.15753e16 1.55161e17

   M pow 0
1 0
0 1

   M&x^:30 =i.2                      Power by repeated multiplication
_1.11401e16 4.15753e16
_4.15753e16 1.55161e17

   x/ 30# ,:M                        Insert x between 30 copies of M
_1.11401e16 4.15753e16
_4.15753e16 1.55161e17

The three techniques for computing matrix powers are applied for various n, with the following timings (milliseconds; J3.01 on Windows 3.1 on 80486/50):

                     30      60      90      12

M&x^:n =i.2        10.4    17.6    26.4    33.5
x/n#,:M             9.3    17.5    25.8    34.0
M pow n             6.6     9.4     9.9     9.3

It remains to define a function to compute the n-th element of the sequence:

   seq=: [: {. M&pow x A0"_

   seq"0 i.10
1 2 7 26 97 362 1351 5042 18817 70226

   seq 30
7.20106e16

The preceding techniques are applicable to any linear recurrence (on any number of terms). For example, the Fibonacci sequence obtains from M=: 0 1,:1 1 and A0=: 0 1.