Расстояние Махаланобиса от инверсии ковариационной матрицы

Я пишу функцию, чтобы взять расстояние Махаланобиса между двумя векторами. Я понимаю, что это достигается с помощью уравнения a ‘* C ^ -1 * b, где a и b – векторы, а C – ковариационная matrix. Мой вопрос в том, есть ли эффективный способ найти обратную матрицу без использования исключения Гаусса-Йордана, или нет никакого способа обойти это? Я ищу способ сделать это сам, а не с любыми предопределенными функциями.

Я знаю, что C – эрмитова положительно определенная matrix, так можно ли каким-либо образом алгоритмически воспользоваться этим фактом? Или есть какой-то умный способ вычислить расстояние Махаланобиса, не вычисляя обратную ковариацию? Любая помощь будет оценена по достоинству.

*** Изменить: неверное уравнение расстояния Махаланобиса указано неверно. Это должно быть x ‘* C ^ -1 * x, где x = (ba), а b и a – два вектора, расстояние которых мы пытаемся найти (спасибо LRPurser). Таким образом, решение, выбранное в выбранном ответе, выглядит следующим образом:

d = x ‘* b, где b = C ^ -1 * x C * b = x, поэтому разрешите для b с использованием факторизации LU или факторизации LDL.

Вы можете (и должны!) Использовать декомпозицию LU вместо явной инвертирования матрицы: решение C x = b с использованием разложения имеет лучшие числовые свойства, чем вычисление C^-1 и умножение вектора b .

Поскольку ваша matrix симметрична, декомпозиция LU эффективно эквивалентна декомпозиции LDL * , которую вы действительно должны использовать в вашем случае. Поскольку ваша matrix также положительно определена, вы должны иметь возможность выполнять это разложение без поворота.


Изменить: обратите внимание, что для этого приложения вам не нужно решать всю проблему C x = b .

Вместо этого, учитывая C = LDL* и разностный вектор v = ab , разрешите L* y = v для y (что в два раза больше, чем полный решатель LU).

Тогда y^t D^-1 y = v^t C^-1 v можно вычислить в линейном времени.

Первое расстояние Махаланобиса (MD) является нормированным расстоянием относительно неопределенности в измерении двух векторов. Когда C=Indentity matrix , MD сводится к евклидову расстоянию, и, следовательно, произведение сводится к векторной норме. Также MD всегда положительно определен или больше нуля для всех ненулевых векторов. По вашей формулировке с соответствующим выбором векторов a и b a*C^-1*b может быть меньше нуля. Надеемся, что разность векторов, которые вы ищете, равна x=(ab) что делает вычисление x^t*C^-1*x где x^t – транспонирование вектора x . Также обратите внимание, что MD=sqrt(x^t*C^-1*x) Так как ваша matrix симметрична и положительно определена, вы можете использовать разложение Холески (MatLab-chol) которое использует половину операций как LU и численно больше стабильный. chol(C)=L где C=L*L^t где L – нижняя треугольная matrix, а L^t – транспонирование L которая делает ее верхней треугольной. Ваш алгоритм должен идти примерно так:

(Matlab)

  x=ab; L=chol(C); z=L\x; MD=z'*z; MD=sqrt(MD);