Матричное умножение с Numpy

Question

Предположим, что у меня есть матрица аффинности A и диагональная матрица D. Как вычислить матрицу Лапласа в Python с помощью nympy?

L = D^(- 1/2) AD^(1/2)

В настоящее время я использую L = D ** (- 1/2) * A * D ** (1/2). Правильно ли это?

Спасибо.

Answer 1

Numpy позволяет экспоненциируется в диагональной матрицы с положительными элементами непосредственно:

m = diag(range(1, 11)) 
print m**0.5 

Однако, это действительно не позволяет экспоненциируются любой матрица непосредственно:

m = matrix([[1, 1], [1, 2]]) 
print m**0.5 

создает наблюдаемый типError (исключение говорит о том, что показатель должен быть целым числом, даже для матрицы, которые можно диагонализировать с положительными коэффициентами).

Итак, до тех пор, пока ваша матрица D диагональна, вы должны иметь возможность напрямую использовать формулу.

Answer 2

Ну, единственная проблема, которую я вижу, заключается в том, что если вы используете Python 2.6.x (без from __future__ import division), то 1/2 будет интерпретироваться как 0, потому что он будет считаться целым делением. Вы можете обойти это, используя вместо этого D ** (-. 5) * A * D **. 5. Вы также можете заставить float-деление с 1./2 вместо 1/2.

Кроме этого, это выглядит правильно для меня.

Edit:

Я пытался экспоненциируются в Numpy массив, а не матрицу до, которая работает с D**.5. Вы можете усилить матричный элемент с помощью numpy.power. Таким образом, вы бы просто использовать

from numpy import power 
power(D, -.5) * A * power(D, .5) 

Answer 3

Имеет ли numpy функцию квадратного корня для матриц? Тогда вы могли бы сделать SQRT (D) вместо (D ** (1/2))

Может быть, формула действительно должна быть написана

L = (D**(-1/2)) * A * (D**(1/2)) 

на основе предыдущего комментария эта формула должна работать в случае D будучи диагональной матрицей (у меня нет шансов доказать это сейчас).

Answer 4

Обратите внимание, что рекомендуется использовать номер array от numpy вместо matrix: см. this paragraph в руководстве пользователя. Путаница в некоторых ответах является примером того, что может пойти не так ... В частности, D ** 0,5 и продукты elementwise, если применяются к массивам numpy, что дало бы вам неправильный ответ. Например:

import numpy as np 
from numpy import dot, diag 
D = diag([1., 2., 3.]) 
print D**(-0.5) 
[[ 1.     Inf   Inf] 
[  Inf 0.70710678   Inf] 
[  Inf   Inf 0.57735027]] 

В вашем случае матрица является диагональной, и поэтому квадратный корень из матрицы просто другая диагональная матрица с квадратным корнем из диагональных элементов.Использование Numpy массивов, уравнение становится

D = np.array([1., 2., 3.]) # note that we define D just by its diagonal elements 
A = np.cov(np.random.randn(3,100)) # a random symmetric positive definite matrix 
L = dot(diag(D**(-0.5)), dot(A, diag(D**0.5)))