Умножение матриц с помощью Numpy

Question

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

L = D ^ (- 1/2) нашей эры ^ (1 /2)

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

Спасибо.

Answer 1

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

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)))

Answer 2

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

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

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

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

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

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

Answer 3

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

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

Редактировать:

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

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

Answer 4

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

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

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

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