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

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

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

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

Спасибо.

4 5
python numpy matrix-multiplication

4 ответа:

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

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

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

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

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

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

Обратите внимание, что вместо matrix рекомендуется использовать array numpy: см. Этот пункт в руководстве пользователя. Путаница в некоторых ответах - пример того, что может пойти не так... В частности, D * * 0.5 и продукты являются элементарно , если их применить к массивам 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)))
4

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

В остальном, по-моему, все правильно.

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

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

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

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

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

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

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

0