Эффективное вычисление обратной матрицы квадратного корня

15

Распространенной проблемой в статистике является вычисление обратного корня квадратного от симметричной положительно определенной матрицы. Что было бы наиболее эффективным способом вычисления этого?

Я натолкнулся на некоторую литературу (которую я еще не читал), и некоторый случайный код R здесь , который я воспроизведу здесь для удобства

# function to compute the inverse square root of a matrix
fnMatSqrtInverse = function(mA) {
  ei = eigen(mA)
  d = ei$values
      d = (d+abs(d))/2
      d2 = 1/sqrt(d)
      d2[d == 0] = 0
      return(ei$vectors %*% diag(d2) %*% t(ei$vectors))
}

Я не совсем уверен, что понимаю линию d = (d+abs(d))/2. Есть ли более эффективный способ вычисления обратного корня квадратного корня матрицы? Функция R eigenвызывает LAPACK .

tchakravarty
источник
d(d+|d|)/2Максимум(d,0)A-1/2A-1/2Икс
@DanielShapero Спасибо за ваш комментарий. Так что, если у меня есть матрица PSD, мне не нужна эта строка? Мое приложение требует вычисления квадратичных форм, таких как . A-1/2ВA-1/2
Чакраварти
Я не знаком с R, но учитывая строку 7, я предполагаю, что она имеет логическое индексирование, как Matlab. Если это так, я предлагаю вам переписать строку 5 как d[d<0] = 0, что более выразительно.
Федерико Полони
Этот код правильный? Я запустил его на простом примере в matlab и нашел ответ неверным. Моя матрица положительно определена, но определенно не симметрична. Пожалуйста, смотрите мой ответ ниже: я передал код в Matlab.
Рони

Ответы:

10

Код , который вы опубликовали ИСПОЛЬЗУЕТ собственное значение разложение симметричной матрицы к вычислительному . A-1/2

Заявление

d = (D + абс (г)) / 2

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

Если вам действительно требуется инверсия квадратного корня симметричной матрицы из , и A достаточно мало (не больше, чем, скажем, 1000 на 1000), то это примерно так же хорошо, как любой другой метод, который вы можете использовать. AA

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

Брайан Борхерс
источник
6
AAAзнак равноВTВВВрA
5

По моему опыту, метод Хайама с полярным Ньютоном работает намного быстрее (см. Главу 6 « Функции матриц » Н. Хайама). В этой короткой заметке есть графики, которые сравнивают этот метод с методами первого порядка. Кроме того, приводятся ссылки на несколько других подходов матрица-квадрат-корень, хотя в основном полярная итерация Ньютона, кажется, работает лучше (и избегает выполнения вычислений по собственным векторам).

% compute the matrix square root; modify to compute inverse root.
function X = PolarIter(M,maxit,scal)
  fprintf('Running Polar Newton Iteration\n');
  skip = floor(maxit/10);
  I = eye(size(M));
  n=size(M,1);
  if scal
    tm = trace(M);
    M  = M / tm;
  else
    tm = 1;
  end
  nm = norm(M,'fro');

  % to compute inv(sqrt(M)) make change here
  R=chol(M+5*eps*I);

  % computes the polar decomposition of R
  U=R; k=0;
  while (k < maxit)
    k=k+1;
    % err(k) = norm((R'*U)^2-M,'fro')/nm;
    %if (mod(k,skip)==0)
    %  fprintf('%d: %E\n', k, out.err(k));
    %end

    iU=U\I;
    mu=sqrt(sqrt(norm(iU,1)/norm(U,1)*norm(iU,inf)/norm(U,inf)));
    U=0.5*(mu*U+iU'/mu);

   if (err(k) < 1e-12), break; end
  end
  X=sqrt(tm)*R'*U;
  X = 0.5*(X+X');
end
suvrit
источник
0

Оптимизируйте свой код:

Вариант 1 - Оптимизируйте свой код R:
a. Вы можете apply()выполнить функцию, dкоторая будет как max(d,0)и d2[d==0]=0в одном цикле.
б. Попробуйте работать ei$valuesнапрямую.

Вариант 2. Использование C ++:
переписать всю функцию в C ++ с помощью RcppArmadillo. Вы все еще сможете позвонить из R.

мощность
источник