Как я могу сделать разложение собственной обобщенной матрицы Фибоначчи до произвольной точности? ⇐ Python

[Anonymous]

Это продолжение моего предыдущего вопроса. Таким образом, я не могу использовать фиксированные библиотеки точности, такие как Numpy. Последовательности Fibonacci, и матрицы для них являются диагонализуемыми, N -й термин может быть эффективно вычислен с использованием матричной диагонализации, и я определил этот метод, чтобы быть близким к оптимальному. означает. Но есть ссылка на итерацию власти, которую я понимаю. Сначала нам нужно умножить матрицу на транспонированную вектор, а затем сам вектор, а затем разделить результат на произведение вектора, транспонированного и самого вектора. Матричный векторный продукт - это вектор длины N, а продукт Vector Dot - это элементы, умножение. Теперь я обнаружил, что транспонирование плоского вектора является самого вектора, и по причинам, которые я не могу объяснить, квадрат вектора собственного вектора всегда 1, поэтому формула упрощается по этому поводу: vec @ mat * vec
Теперь нам нужно найти норму и переработать новое вектор, а затем и приносят на квадраты. Python: sum (i * i для i in vec) ** 0.5 .
Затем мы разделяем вектор на норму и повторяем весь процесс.
и инициализируем вектор, чтобы быть всеми, поэтому это гарантированно неортогонально для eigen. />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
< /code>
Теперь, используя отношения, которые я уже обнаружил, столбцы в следующей итерации состояния матрицы Фибоначчи сдвинуты на 1, а первый элемент заменяется суммой вектора, поэтому полностью оптимизированный код: < /p>
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
< /code>
Но я хорошо знаю, что мой код не делает всю матричную разложение, используя Numpy в качестве ссылки, я обнаружил: < /p>
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)
< /code>
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 ]]))
< /code>
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?

Подробнее здесь: https://stackoverflow.com/questions/795 ... rbitrary-p