Быстрее мощность матрицы, чем numpy?

Question

Мне нужно было вычислить Q^N для множества различных значений N (от 1 до 10000) и Numpy было слишком медленным.

Я спросил у math.stackexchange.com, если бы я мог избежать вычисления Q^N для моей конкретной потребности, и кто-то ответил мне, что вычисление Q^N должно быть довольно быстрым, используя метод P D^N P^-1.

Так в основном, вместо того, чтобы делать:

import numpy as np 
from numpy import linalg as LA 
... 
LA.matrix_power(m, N) 

Я пробовал:

diag, P = LA.eig(m) 
DN = np.diag(diag**N) 
P1 = LA.inv(P) 

P*DN*P1 

И я получаю ту же матрицу, как результат (пытался на одном примере)

О более сложной матрице, Q:

% timeit.Timer('Q**10000', setup=setup).repeat(2, 100) 
[5.87254786491394, 5.863131046295166] 

% timeit.Timer('diag, P = linalg.eig(Q); DN=np.diag(diag**10000);P1=linalg.inv(P); P*DN*P1', setup=setup).repeat(2, 100) 
[2.0032401084899902, 2.018735885620117] 

И что касается моей первоначальной проблемы, второй метод позволяет мне вычислить P, diag and P1 только один раз и использовать его тысячи раз. Этот метод используется в 8 раз быстрее.

Мои вопросы:

• В этом случае не представляется возможным использовать этот последний метод для вычисления Q^N?
• Можно ли использовать его в моем случае (матрица Q, как указано here)?
• Есть ли в numpy функция, которая уже делает это?

Edit:

• Оказывается, что для другой матрицы P не является обратимым. Поэтому я добавляю новый вопрос: как я могу изменить матрицу P, чтобы она стала обратимой, но результирующая матрица не слишком изменена? Я имею в виду, что это нормально, если значения близки к реальному результату, закрываясь, я имею в виду ~ 0,0001.

Answer 1

Я частично отвечая на мой вопрос:

Согласно source code, я думаю, что Numpy использует алгоритмы быстрого возведения в степень:

# binary decomposition to reduce the number of Matrix 
# multiplications for n > 3. 
beta = binary_repr(n) 
Z, q, t = M, 0, len(beta) 
while beta[t-q-1] == '0': 
    Z = N.dot(Z, Z) 
    q += 1 
result = Z 
for k in range(q+1, t): 
    Z = N.dot(Z, Z) 
    if beta[t-k-1] == '1': 
     result = N.dot(result, Z) 
return result 

Который медленнее в моем случае, когда n является большим, чем вычисления собственных значений и собственных векторов и вычислить M^N равным P D^N P^-1.

Теперь, что касается моих вопросов:

В этом случае не представляется возможным использовать этот последний метод для вычисления Q^N?

Когда некоторые собственные значения равны, то не будет возможно инвертировать P. кого предложил сделать это в Numpy на issue tracker. Ответ был: «Ваш подход действителен только для не дефектных плотных матриц».

Можно ли использовать его в моем случае (матрица Q, как указано здесь)?

Не всегда, я мог бы иметь несколько равных собственных значений.

Есть ли в numpy функция, которая уже делает это?

Я думаю, что это в SciPy: https://github.com/scipy/scipy/blob/v0.12.0/scipy/linalg/matfuncs.py#L57

Таким образом, мы также можем сделать это:

LA.expm(n*LA.logm(m)) 

вычислить т^п.

Как изменить матрицу P, чтобы она стала обратимой, но результирующая матрица не слишком изменена? Я имею в виду, что это нормально, если значения близки к реальному результату, закрываясь, я имею в виду ~ 0,0001.

Я не могу просто добавить значение epsilon, потому что метод декомпозиции является разумным, когда значения слишком близки. Я уверен, что это может привести к непредсказуемым ошибкам.

Answer 2

Вы уже выяснили, что ваши собственные значения будут (0, a, b, c, ..., 1). Позвольте мне переименовать ваши параметры, чтобы собственные значения были (0, e1, e2, e3, ..., 1). Для того, чтобы выяснить, собственный вектор (v0, v1, v2, ..., v(n-1)), соответствующий собственному значению ej, вы должны решить систему уравнений:

v1     = v0*ej 
v1*e1 + v2*(1-e1)  = v1*ej 
v2*e2 + v3*(1-e2)  = v2*ej 
... 
vj*ej + v(j+1)*(1-ej) = vj*ej 
... 
v(n-1)    = v(n-1)*ej 

Это более или менее ясно, что, если все ваши ei различны, и никто не равен 0 или 1, то решение хорошо определено всегда, а при работе с ej результирующий собственный вектор имеет первые j компоненты, отличные от нуля, а остальные равны нулю. Это гарантирует, что ни один из собственных векторов не является линейной комбинацией остальных, и, следовательно, матрица собственных векторов обратима.

Проблема возникает, когда некоторые из ваших ei либо 0, либо 1, либо повторяется.Я не был в состоянии придумать доказательство этого, но поэкспериментировать со следующим кодом, кажется, что вы должны беспокоиться только если два из ваших ei одинаковы и отличаются от 1:

>>> def make_mat(values): 
...  n = len(values) + 2 
...  main_diag = np.concatenate(([0], values, [1])) 
...  up_diag = 1 - np.concatenate(([0], values)) 
...  return np.diag(main_diag) + np.diag(up_diag, k=1) 
>>> make_mat([4,5,6]) 
array([[ 0, 1, 0, 0, 0], 
     [ 0, 4, -3, 0, 0], 
     [ 0, 0, 5, -4, 0], 
     [ 0, 0, 0, 6, -5], 
     [ 0, 0, 0, 0, 1]]) 
>>> a, b = np.linalg.eig(make_mat([4,5,6])) 
>>> a 
array([ 0., 4., 5., 6., 1.]) 
>>> b 
array([[ 1.  , 0.24253563, -0.18641093, 0.13608276, 0.4472136 ], 
     [ 0.  , 0.9701425 , -0.93205465, 0.81649658, 0.4472136 ], 
     [ 0.  , 0.  , 0.31068488, -0.54433105, 0.4472136 ], 
     [ 0.  , 0.  , 0.  , 0.13608276, 0.4472136 ], 
     [ 0.  , 0.  , 0.  , 0.  , 0.4472136 ]]) 

И теперь для некоторых тестов:

>>> a, b = np.linalg.eig(make_mat([1,0,3])) # having a 0 or 1 is OK 
>>> b 
array([[ 1.  , 0.70710678, 0.  , 0.  , 0.  ], 
     [ 0.  , 0.70710678, 0.  , 0.  , 0.  ], 
     [ 0.  , 0.  , 1.  , 0.31622777, 0.57735027], 
     [ 0.  , 0.  , 0.  , 0.9486833 , 0.57735027], 
     [ 0.  , 0.  , 0.  , 0.  , 0.57735027]]) 
>>> a, b = np.linalg.eig(make_mat([1,1,3])) # repeating 1 is OK 
>>> b 
array([[ 1.  , 0.70710678, 0.  , 0.  , 0.  ], 
     [ 0.  , 0.70710678, 0.  , 0.  , 0.  ], 
     [ 0.  , 0.  , 1.  , 0.  , 0.  ], 
     [ 0.  , 0.  , 0.  , 1.  , 0.70710678], 
     [ 0.  , 0.  , 0.  , 0.  , 0.70710678]]) 
>>> a, b = np.linalg.eig(make_mat([0,0,3])) # repeating 0 is not OK 
>>> np.round(b, 3) 
array([[ 1. , -1. , 1. , 0.035, 0.447], 
     [ 0. , 0. , 0. , 0.105, 0.447], 
     [ 0. , 0. , 0. , 0.314, 0.447], 
     [ 0. , 0. , 0. , 0.943, 0.447], 
     [ 0. , 0. , 0. , 0. , 0.447]]) 
>>> a, b = np.linalg.eig(make_mat([2,3,3])) # repeating other values are not OK 
>>> np.round(b, 3) 
array([[ 1. , 0.447, -0.229, -0.229, 0.447], 
     [ 0. , 0.894, -0.688, -0.688, 0.447], 
     [ 0. , 0. , 0.688, 0.688, 0.447], 
     [ 0. , 0. , 0. , 0. , 0.447], 
     [ 0. , 0. , 0. , 0. , 0.447]])