Нахождение собственных векторов ковариационной матрицы для создания трехмерной ограничивающей сферы

Question

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

Вот то, с чем я работаю: A) Точки в трехмерном пространстве B) Ковариационная матрица 3x3, хранящаяся в матричном классе 4x4 (на который ссылаются ячейки m0,m1,m2,m3,m4,ect; вместо строк и столбцов)

Я нашел 3 собственных значения для ковариационной матрицы точек и настроил функцию для преобразования матрицы в уменьшенную форму ряда строк (rref) посредством исключения Гаусса.

Я протестировал обе эти функции на рисунках в примерах, которые я нашел в Интернете, и они, кажется, работают правильно.

Следующим шагом является поиск собственных векторов с использованием уравнения:(M - λ * I) * V

... где M - ковариационная матрица, λ - одно из собственных значений, I - единичная матрица, а V - собственный вектор.

Тем не менее, я, похоже, неправильно строю матрицу 4x3 до ее наведения, поскольку крайний правый столбец, где должны быть рассчитаны компоненты собственного вектора, равен 0 до и после выполнения rref. Я понимаю, почему они равны нулю после (без каких-либо констант, самое простое решение линейной системы уравнений - это все коэффициенты, равные нулю), но я не знаю, что там поставить.

Вот функция до сих пор:

Vect eigenVector (const Matrix & M, const float eval) {
   Матрица А = Матрица (М);
   A -= Матрица (IDENTITY)*eval;
   A.rref();
   возврат Vect(A[м3],A[м7],A[м11]);
}

Ковариационная матрица 3x3 передается как M, а собственное значение как eval. Matrix(IDENTITY) возвращает единичную матрицу. m3,m7 и m11 соответствуют крайнему правому столбцу матрицы 4x3.

Вот пример матрицы 3х3 (хранится в классе матрицы 4х4), которую я использую для проверки функций:

Матрица (1,5f, 0,5f, 0,75f, 0,
       0,5f, 0,5f, 0,25f, 0,
      0,75f, 0,25f, 0,5f, 0,
          0,     0,    0, 0);

Я правильно (?) Получаю собственные значения 2.097, 0.3055, 0.09756 из моей другой функции.

eigenVector () выше корректно вычитает переданное собственное значение из диагонали (0,0 1,1 2,2)

Матрица А после rref ():

[(1, 0, 0, -0),
(-0, 1, 0, -0),
(-0, -0, 1, -0),
(0, 0, 0, -2.09694)]

Для функции rref () я использую переведенную функцию python, найденную здесь: http://elonen.iki.fi/code/misc-notes/python-gaussj/index.html

Как должна выглядеть матрица, которую я передаю rref (), чтобы получить собственный вектор?

Спасибо

Answer 1

(M - λI)V - это не уравнение, это просто выражение. Однако (M - λI)V = 0 есть. И это уравнение, которое связывает собственные векторы с собственными значениями.

Итак, если ваш rref функция работает, я хотел бы представить, что вы создаете расширенную матрицу как [(M - λI) | 0], где 0 обозначает нулевой вектор. Это похоже на то, что вы уже делаете, поэтому я должен предположить, что ваш rref функция нарушена. Или, в качестве альтернативы, он не знает, как обращаться с матрицами 4х4 (в отличие от матриц 4х3, чего и следует ожидать от расширенной матрицы).

Answer 2

Ах, после нескольких часов изнурительных исследований, мне удалось решить мою проблему.

Проблема в том, что не существует "одного" набора собственных векторов, а существует бесконечное число с различными величинами.

Метод, который я выбрал, состоял в том, чтобы использовать REF (форму эшелона строки) вместо RREF, оставляя достаточно информации в матрице, чтобы позволить мне заменить произвольное значение на z, и работать в обратном направлении, чтобы найти y и x. Затем я нормализовал вектор, чтобы получить собственный вектор единицы, который должен работать для моих целей.

Мой окончательный код:

Vect eigenVector (const Matrix & M, const float eVal) {
   Матрица А = Матрица (М);
   A -= Матрица (IDENTITY)*eVal;
   A.ref();
   поплавок К = 16; // Произвольное значение
   поплавок J = -K*A[м6]; // подставим в K, чтобы найти J
   поплавок I = -K*A[м2]-J*A[м1]; // Подставим в K и J, чтобы найти I

   Vect eVec = Vect(I,J,K);
   eVec.norm(); // Нормализуем собственный вектор

   вернуть eVec;
}

Единственная странность заключается в том, что собственные векторы выходят лицом в противоположном направлении, чем я ожидал (они были сведены на нет!), Но это спорная проблема.