[Python] How do I do Eigen decomposition of a generalized Fibonacci matrix to arbitrary...

This is a follow up to my previous question.

I want to efficiently compute Nth term of higher order generalized Fibonacci numbers, N is sufficiently large such that the Nth term is guaranteed to be outside of UINT64 range. So I cannot use fixed precision libraries like NumPy.

It is well known equations with exponents larger than 4 (quintic equations and above) can't be solved algebraically, so we cannot leverage closed form characteristic polynomials of generalized Fibonacci numbers of order 5 and above.

But matrix exponentiation is guaranteed to work for higher order Fibonacci sequences, and the matrices for them are diagonalizable, the Nth term can be computed efficiently using matrix diagonalization, and I have determined this method to be close to optimal.

So I need to do Eigen decomposition of a matrix numerically to arbitrary precision.

I have no idea what the notation used means. But there is a link to Power iteration which I understand.

Basically, start with a vector of length N that isn't orthogonal to the Eigen vector of the square matrix (whose side is N). First we need to multiply the matrix by the vector transposed and then the vector itself, then divide the result by the product of the vector transposed and the vector itself. Matrix vector product is a vector of length N, and vector dot product is itemwise multiplication.

So in Python syntax, this is: vec.T @ mat * vec / (vec.T * vec). Now I discovered that the transposition of a flat vector is the vector itself, and for reasons I cannot explain, the square of a eigen vector is always 1, so the formula is simplified to this: vec @ mat * vec

Now we need to find the norm and renormalize the new vector, it is such a fancy name for summing all the squares of the numbers in the vector and then taking the square root, or in Python: sum(i * i for i in vec) ** 0.5.

Then we divide the vector by the norm and repeat the whole process.

And I initialize the vector to be all ones so it is guaranteed to be non-orthogonal to the eigen vector.

So, here is the code:

import gmpy2
import numpy as np

Matrix = list[int]

ONACCI_MATRICES = {}

def onacci_matrix(n: int) -> Matrix:
if matrix := ONACCI_MATRICES.get(n):
return matrix

mat = [1] * n + [0] * (n * (n - 1))
for i in range(1, n):
mat[i * n + i - 1] = 1

ONACCI_MATRICES[n] = mat
return mat


def mat_vec_mult(mat: Matrix, vec: Matrix, side: int) -> Matrix:
return [
sum(a * b for a, b in zip(mat[i : i + side], vec))
for i in range(0, side * side, side)
]


gmpy2.get_context().precision = int(100 / gmpy2.log10(2) + 1)


def power_iteration(mat: Matrix, k: int, lim: int) -> Matrix:
vec = [gmpy2.mpfr(1)] * k
for _ in range(lim):
vec1 = mat_vec_mult(mat, vec, k)
d = gmpy2.sqrt(sum(i * i for i in vec1))
vec = [i / d for i in vec1]

vec1 = mat_vec_mult(mat, vec, k)
return sum(a * b for a, b in zip(vec, vec1)), vec


Now, using the relationships I have already discovered, the columns in the next iteration of Fibonacci matrix state are shifted down by 1, and the first item is replaced by the sum of the vector, so the fully optimized code is:

def onacci_eigen(k: int, lim: int) -> Matrix:
vec = [gmpy2.mpfr(1)] * k
for _ in range(lim):
vec1 = [sum(vec)] + vec[:-1]
norm = gmpy2.sqrt(sum(i * i for i in vec1))
vec = [i / norm for i in vec1]

vec1 = [sum(vec)] + vec[:-1]
return sum(a * b for a, b in zip(vec, vec1)), vec


But I am well aware of the fact my code doesn't do the whole matrix decomposition, using NumPy as reference, I have found:

def onacci_eigen_np(k: int):
mat = np.zeros((k, k), dtype=np.uint8)
mat[0] = 1
for i in range(1, k):
mat[i, i - 1] = 1

return np.linalg.eig(mat)


In [177]: onacci_eigen(2, 16)
Out[177]:
(mpfr('1.61803398874988959589659781966119622236443053427999997832557479220999483606819424403126969324727864072',333),
[mpfr('0.850650782872223955022787840306148985083417492530648086674174311763827379148760465636071429885324565541',333),
mpfr('0.525731153346339650221127005018931861137081167016813078335393334321529537345422006045203588826185499676',333)])

In [178]: onacci_eigen_np(2)
Out[178]:
EigResult(eigenvalues=array([ 1.61803399, -0.61803399]), eigenvectors=array([[ 0.85065081, -0.52573111],
[ 0.52573111, 0.85065081]]))

In [179]: onacci_eigen(3, 16)
Out[179]:
(mpfr('1.839286653690119133028569875712653669087335320613832305625305015512987239734543867343593384912450164',333),
[mpfr('0.850340279862825348236116977938202426830832418712093774829379207962850131102785349206559081374605648379',333),
mpfr('0.462320602298743499525549933221019435813052948664915507101373253172318146117945675600021936374904261574',333),
mpfr('0.251358447506621329352814452598084994094648576600430235326923413674553183407741581212001741035260789156',333)])

In [180]: onacci_eigen_np(3)
Out[180]:
EigResult(eigenvalues=array([ 1.83928676+0.j , -0.41964338+0.60629073j,
-0.41964338-0.60629073j]), eigenvectors=array([[-0.85034021+0.j , -0.14119411-0.37520324j,
-0.14119411+0.37520324j],
[-0.46232063+0.j , -0.30942518+0.44705011j,
-0.30942518-0.44705011j],
[-0.25135865+0.j , 0.73735271+0.j ,
0.73735271-0.j ]]))

In [181]: onacci_eigen(4, 16)
Out[181]:
(mpfr('1.92756223342068232349610814445895076913312005500745526311521135939897015221510102277085793961458193032',333),
[mpfr('0.857153590513098422301230530952927225060858068163042414703783315715733220474264891990176420761499548907',333),
mpfr('0.444682407921392626174463658040140355928890877348872363567201899637692039777354781848554832843424133843',333),
mpfr('0.230696795811217076801604375082105672439567442961537961628103903254219452734144738765055900155347524604',333),
mpfr('0.119684028834984385725222174819194615542649649688339912067213468714238731999300222194281164844625699865',333)])

In [182]: onacci_eigen_np(4)
Out[182]:
EigResult(eigenvalues=array([ 1.92756198+0.j , -0.07637893+0.81470365j,
-0.07637893-0.81470365j, -0.77480411+0.j ]), eigenvectors=array([[ 0.85715349+0.j , -0.09751756+0.33857661j,
-0.09751756-0.33857661j, 0.31523984+0.j ],
[ 0.44468271+0.j , 0.42308563+0.08003246j,
0.42308563-0.08003246j, -0.40686392+0.j ],
[ 0.23069697+0.j , 0.04911753-0.5239171j ,
0.04911753+0.5239171j , 0.52511843+0.j ],
[ 0.1196833 +0.j , -0.6430769 +0.j ,
-0.6430769 -0.j , -0.67774348+0.j ]]))

In [183]: onacci_eigen(5, 16)
Out[183]:
(mpfr('1.9659474114747634239368013905053543613908549162474134872580997206921779654734866152502941435672985726',333),
[mpfr('0.861467215187968428723649888178636920487027712949253173373040938936655670698160778471770465480475619984',333),
mpfr('0.438193862570299377379587205660059318052136742144365550315894896468306717641558127893410820652331043709',333),
mpfr('0.222891136965451570141305174696044140066685007875584632245869788254444544661357278766279850404309764097',333),
mpfr('0.113375681300648031827460418686209579835989861656615983444017189920553612517748621742739205364731840502',333),
mpfr('0.057670372929869948367396972579646922674674571679137745911613868271570224905410941730481509887862842752',333)])

In [184]: onacci_eigen_np(5)
Out[184]:
EigResult(eigenvalues=array([ 1.96594824+0.j , 0.19537659+0.84885364j,
0.19537659-0.84885364j, -0.67835071+0.45853619j,
-0.67835071-0.45853619j]), eigenvectors=array([[-0.86146684+0.j , -0.20189765+0.25699988j,
-0.20189765-0.25699988j, 0.20038647+0.19197533j,
0.20038647-0.19197533j],
[-0.43819406+0.j , 0.23553899+0.29206031j,
0.23553899-0.29206031j, -0.07145523-0.33130376j,
-0.07145523+0.33130376j],
[-0.22289196+0.j , 0.3874071 -0.18831127j,
0.3874071 +0.18831127j, -0.15429702+0.38409776j,
-0.15429702-0.38409776j],
[-0.11337631+0.j , -0.11092093-0.48191871j,
-0.11092093+0.48191871j, 0.41883037-0.28311149j,
0.41883037+0.28311149j],
[-0.05767004+0.j , -0.56772886+0.j ,
-0.56772886-0.j , -0.61742453+0.j ,
-0.61742453-0.j ]]))

In [185]: onacci_eigen(6, 16)
Out[185]:
(mpfr('1.98358305374900825496807475402132087846852695358925128301033942225440908180872792292068910732815562937',333),
[mpfr('0.863738989414237627573489167420515042155134684491856711393028223100192381233871590661405626679603882099',333),
mpfr('0.435444052359025770906718098292102305385375052195794885332245444423951113107552568689013662941192744191',333),
mpfr('0.219524895859618218306890496627457403937545340156672314076372625231620069403920826656866625559799029342',333),
mpfr('0.110672386872126423376689893976513932549903736608000832702793574988932501093306951193633774890959346395',333),
mpfr('0.0557925802575993100993804234733299323919677614855456275102808038249108520700391223142538793454193395756',333),
mpfr('0.0281259008899451455546211036328818000809421872502582926611627952215610951244247623006821964670892534942',333)])

In [186]: onacci_eigen_np(6)
Out[186]:
EigResult(eigenvalues=array([ 1.98358284+0.j , 0.39029203+0.81786166j,
0.39029203-0.81786166j, -0.8403091 +0.j ,
-0.46192891+0.71914438j, -0.46192891-0.71914438j]), eigenvectors=array([[ 0.86373937+0.j , -0.24604208+0.18916884j,
-0.24604208-0.18916884j, -0.24267103+0.j ,
0.0727323 +0.24663839j, 0.0727323 -0.24663839j],
[ 0.43544406+0.j , 0.07146109+0.33493779j,
0.07146109-0.33493779j, 0.28878781+0.j ,
0.19679976-0.2275479j , 0.19679976+0.2275479j ],
[ 0.21952401+0.j , 0.3675281 +0.08801268j,
0.3675281 -0.08801268j, -0.34366856+0.j ,
-0.34843385-0.04984825j, -0.34843385+0.04984825j],
[ 0.11067045+0.j , 0.26232195-0.32419412j,
0.26232195+0.32419412j, 0.40897874+0.j ,
0.17124646+0.37451473j, 0.17124646-0.37451473j],
[ 0.05579321+0.j , -0.19819618-0.41532249j,
-0.19819618+0.41532249j, -0.48670036+0.j ,
0.26038911-0.40538135j, 0.26038911+0.40538135j],
[ 0.02812749+0.j , -0.50781509+0.j ,
-0.50781509-0.j , 0.57919206+0.j ,
-0.56369954+0.j , -0.56369954-0.j ]]))

In [187]: onacci_eigen(7, 16)
Out[187]:
(mpfr('1.99196215673683730102465372048514883373603183657729478919904292662019648954693003866195864703727267182',333),
[mpfr('0.864885515585064461373786114718374480355277984736331348986589266091967623451587299263112217554012928619',333),
mpfr('0.43418873989287852678012919397803697628504247534472892445357785344244700208211704404295277823149688493',333),
mpfr('0.217968496225546316700722480731875638209296713544896128331742164500987256941168413891941561729505024996',333),
mpfr('0.10942244648126059026972999767079867615110982909296393834315932872190766714258734754597020505814217594',333),
mpfr('0.054931457653807916358589192921262048611859385090989947385035415178193013629243438737751652512277291872',333),
mpfr('0.0275769812624266740551965089369946615871556622019990039240572105792482575655447965160259391342531200748',333),
mpfr('0.0138452520779086002295890261401114314472645503135015582026099843662321648801821314963548638926951692937',333)])

In [188]: onacci_eigen_np(7)
Out[188]:
EigResult(eigenvalues=array([ 1.9919642 +0.j , 0.52886126+0.76534196j,
0.52886126-0.76534196j, -0.24065634+0.849197j ,
-0.24065634-0.849197j , -0.78418702+0.36004972j,
-0.78418702-0.36004972j]), eigenvectors=array([[ 0.86488565+0.j , 0.26349599-0.13933657j,
0.26349599+0.13933657j, 0.02103513-0.24366028j,
0.02103513+0.24366028j, -0.189281 -0.11841467j,
-0.189281 +0.11841467j],
[ 0.43418734+0.j , 0.03779881-0.31816583j,
0.03779881+0.31816583j, -0.27209748+0.05233986j,
-0.27209748-0.05233986j, 0.14208791+0.216241j ,
0.14208791-0.216241j ],
[ 0.21796945+0.j , -0.2582671 -0.22785406j,
-0.2582671 +0.22785406j, 0.14110608+0.28042893j,
0.14110608-0.28042893j, -0.04508011-0.29644979j,
-0.04508011+0.29644979j],
[ 0.10942438+0.j , -0.35932337+0.08915608j,
-0.35932337-0.08915608j, 0.26208981-0.2404385j ,
0.26208981+0.2404385j , -0.09587276+0.33401577j,
-0.09587276-0.33401577j],
[ 0.05493291+0.j , -0.14073415+0.37224475j,
-0.14073415-0.37224475j, -0.34304972-0.21141469j,
-0.34304972+0.21141469j, 0.26248775-0.30542094j,
0.26248775+0.30542094j],
[ 0.02757726+0.j , 0.24318852+0.35193044j,
0.24318852-0.35193044j, -0.12447918+0.43924604j,
-0.12447918-0.43924604j, -0.42413717+0.19473731j,
-0.42413717-0.19473731j],
[ 0.01384425+0.j , 0.45983425+0.j ,
0.45983425-0.j , 0.5172487 +0.j ,
0.5172487 -0.j , 0.54086226+0.j ,
0.54086226-0.j ]]))


So my code finds the first Eigen value and the first column of Eigen vectors, and for some odd powers the first column of the Eigen vectors have negative signs but not all (7 is a counter example), which is odd.

How do I do the full Eigen decomposition of a matrix using my existing code?

Continue reading...