GSL / BLAS: умножить матрицу на обратную матрицу - PullRequest
Новая
       21

GSL / BLAS: умножить матрицу на обратную матрицу

4 голосов
/ 23 августа 2010

Я использую GNU GSL для выполнения некоторых матричных вычислений. Я пытаюсь умножить матрицу B на обратную матрицу A.

Теперь я заметил, что BLAS-часть GSL имеет функцию для этого, но только если A треугольная. Есть ли конкретная причина для этого? Кроме того, что было бы самым быстрым способом сделать это вычисление? Должен ли я инвертировать A, используя разложение LU, или есть лучший способ?

FWIW, A имеет вид P ' G P, где P - нормальная матрица, P' - ее обратная, а G - диагональная матрица.

Спасибо большое :))

Ответы [ 2 ]

4 голосов
/ 05 сентября 2010

Я считаю, что Адриен прав, потому что у BLAS нет обратной функции для квадратных матриц.Это зависит от матрицы, которую вы используете для оптимизации исчисления ее обратного.

В общем случае вы можете использовать разложение LU, которое справедливо для любой квадратной матрицы.То есть, что-то вроде: 

gsl_linalg_LU_decomp(A, p, signum);
gsl_linalg_LU_invert(A, p, invA);

, где A - квадратная матрица, которую вы хотите инвертировать, p - это gsl_permutation объект (объект перестановки, где закодирована матрица перестановок), signum -знак перестановки и invA является обратным к A.

Поскольку вы утверждаете, что A = P' G P является P нормальным и G диагональю, вероятно, A является нормальной матрицей.Я не использовал их некоторое время, но для них должна быть теорема о факторизации, которую вы можете найти в главе 14 справочного руководства GSL или даже лучше, в книге по линейной алгебре.

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

Надеюсь, это поможет!

3 голосов
/ 23 августа 2010

в двух словах:

Тот факт, что поддерживаются только треугольные матрицы, объясняется просто тем, что эту операцию действительно легко выполнить для треугольной матрицы (ее называют обратной подстановкой), а BLAS предоставляет подпрограммы только для функций низкого уровня. Любая функция более высокого уровня обычно находится в LAPACK, который использует BLAS для блочных операций.

В конкретном случае, с которым вы имеете дело: A:=P'*G*P, если P - нормальная матрица, то обратное значение A может быть записано inv(A) := P'*inv(G)*P. Затем у вас есть два варианта:

  1. Построить inv(A) и умножить его на B.
  2. последовательно умножить B на различные коэффициенты inv(A): умножить B на P (слева), затем изменить масштаб каждой строки результата с инверсией диагональных элементов G и затем умножить снова с P' (снова слева).

В зависимости от конкретного случая вы можете выбрать свой метод. В любом случае вам понадобятся dgemm (умножение матрицы на матрицу, BLAS уровня 3) и dscal (масштабирование линий, BLAS уровня 1), при условии двойной точности.

Надеюсь, это поможет.

A.

...