Essays/Linear Recurrences

From J Wiki
Jump to: navigation, search

Linear Recurrences and Matrix Powers

I read with enjoyment and admiration accounts of Eugene McDonnell's excellent adventure in Heron's Rule & Integer-Area Triangles in the January Vector. The article discusses Sloane sequence A001075, 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,

   mp=: +/ .*   NB. Matrix multiplication
   M=: 0 1,:_1 4
   A0=: 1 2
   M mp A0
2 7
   M mp M mp A0                    (M mp M) mp A0
7 26				7 26
   M mp M mp M mp A0               (M mp M mp M) mp A0
26 97				26 97
   M&mp^:(i.5) A0                  (M&mp^:(i.5) =i.2) mp 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                   NB. Binary representation of 30
1 1 1 1 0
   I.|.b                        NB. Powers corresponding to 1s in b
1 2 3 4

   M mp M			NB. Squaring a matrix
_1  4
_4 15
   mp~ M
_1  4
_4 15
   mp~^:1 2 3 4 M               NB. Powers with squaring as the function
        _1          4
        _4         15

       _15         56
       _56        209

     _2911      10864
    _10864      40545

_109552575  408855776
_408855776 1525870529

   (M&mp^:(2^1 2 3 4) =i.2) -: mp~^:1 2 3 4 M
1
   mp/ mp~^:1 2 3 4 M           NB. Product of selected powers
_1.11401e16 4.15753e16
_4.15753e16 1.55161e17

   pow=: 4 : 'mp/ mp~^:(I.|.#:y) x'

   M pow 30
_1.11401e16 4.15753e16
_4.15753e16 1.55161e17
   M pow 0
1 0
0 1
   M&mp^:30 =i.2		NB. Power by repeated multiplication
_1.11401e16 4.15753e16
_4.15753e16 1.55161e17
   mp / 30# ,:M			NB. Insert mp 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 (microseconds; J6.01 on an AMD 64 Athlon 3200+ 2.2 GHz):

30 60 90 120
M&mp^:n =i.2 61.9 93.2 124.9 155.3
mp/n#,:M 55.6 89.1 129.5 163.2
M pow n 48.6 67.2 68.3 67.7

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

   seq=: [: {. M&pow mp 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 .



See also



Contributed by Roger Hui. Substantially the same text previously appeared in Vector, Volume 12, Number 4, April 1996.