Essays/Linear Recurrences
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.