Если вы ЗНАЕТЕ, что матрица имеет ранг 3, то ровно 3 преобразования Хаусхолдера будет достаточно, чтобы уменьшить матрицу nxm до матрицы 3xm.Теперь можно легко преобразовать это в проблему собственных значений.Вычислить характеристический полином.Например, рассмотрим эту матрицу (я буду использовать MATLAB для выполнения работы):
>> A = rand(10,3);
>> B = A*rand(3,6);
Очевидно, что B будет иметь ранг 3, но если вы мне не верите, ранг подтверждает это утверждение.
>> rank(B)
ans =
3
>> size(B)
ans =
10 6
Итак, B - матрица 10x6 с рангом 3. SVD также подтверждает этот факт.
>> svd(B)
ans =
6.95958965358531
1.05904552889497
0.505730235622534
7.37626877572817e-16
2.36322066654691e-16
1.06396598411356e-16
Мне лень писать преобразования Домохозяина.(У меня есть некоторый код для этого, но ...) Так что я буду использовать QR, чтобы выручить меня.
>> [Q,R] = qr(B,0);
>> C = Q(:,1:3)'*B
C =
-2.0815 -1.7098 -3.7897 -1.6186 -3.6038 -3.0809
0.0000 0.91329 0.78347 0.44597 -0.072369 0.54196
0.0000 0.0000 -0.2285 -0.43721 -0.85949 -0.41072
Умножение здесь показывает, что мы увидели бы после трех преобразований Домохозяев,Посмотрите, что это верхняя треугольная, как мы и ожидали.Далее вычисляем характеристический полином.Я сделаю это так символически здесь, используя мои собственные инструменты, но вычисления - это всего лишь алгебра.
>> sympoly lambda
>> det(C*C' - lambda*eye(3))
ans =
13.8942 + 66.9996*lambda + 49.8132*lambda^2 + lambda^3
>> P = det(C*C' - lambda*eye(3))
P =
13.8942 - 66.9996*lambda + 49.8132*lambda^2 - lambda^3
Каковы корни P, так что собственные значения C * C '?
>> r = roots(P)
r =
48.436
1.1216
0.25576
Мы знаем, что собственные значения должны быть здесь квадратами сингулярных значений, поэтому хороший тест здесь - посмотреть, восстановили ли мы сингулярные значения, найденные svd.Так что, снова расширив формат отображения, мы видим, что это очень хорошо.
>> sqrt(roots(P))
ans =
6.95958965358531
1.05904552889497
0.505730235622533
Математика может быть интересной, когда она работает.Что мы можем сделать с этой информацией?Если я знаю конкретное собственное значение, я могу вычислить соответствующий собственный вектор.По сути, нам нужно решить линейную однородную систему уравнений 3x3
(C*C' - eye(3)*r(3)) * X = 0
Опять же, я буду ленив здесь, чтобы найти решение, фактически не написав никакого кода.Гауссовское исключение сделало бы это, хотя.
>> V = null((C*C' - eye(3)*r(3)))
V =
-0.171504758161731
-0.389921448437349
0.904736084157367
Итак, у нас есть V, собственный вектор C * C '.Мы можем убедить себя, что это один из следующих вызовов svd.
>> svd(C - V*(V'*C))
ans =
6.9595896535853
1.05904552889497
2.88098729108798e-16
Таким образом, вычитая этот компонент C в направлении V, мы получаем матрицу ранга 2.
Точно так же мы можем преобразовать V в исходное проблемное пространство и использовать его для преобразования матрицы B, нашей исходной матрицы, с обновлением ранга один B.
>> U = Q(:,1:3)*V;
>> D = B - U*(U'*B);
Что такое ранг D?
>> svd(D)
ans =
6.95958965358531
1.05904552889497
2.62044567948618e-15
3.18063391331806e-16
2.16520878207897e-16
1.56387805987859e-16
>> rank(D)
ans =
2
Как видите, хотя я много занимался математикой, вызывая svd, QR, rank и т. Д., Все несколько раз, в конце концов, фактические вычисления были относительно тривиальными.Мне нужны были только ...
- 3 домохозяина.(Сохраните их для последующего использования. Обратите внимание, что преобразование Домохозяина - это просто обновление первого уровня вашей матрицы.)
- Вычисление характеристического полинома.
- Вычисление наименьшего корня с использованием вашего любимого метода длякубический полином.
- Восстановить соответствующий собственный вектор.Здесь достаточно исключения Гаусса.
- Верните этот собственный вектор в исходное пространство, используя преобразования домохозяев, которые мы изначально делали.
- Обновление матрицы первого ранга.
Все эти вычислительные шаги будут быстрыми и эффективными для ЛЮБОЙ матрицы размера, если мы знаем, что фактический ранг равен 3. Даже не стоит бумаги по этому вопросу.
РЕДАКТИРОВАТЬ:
Свопрос был пересмотрен так, что сама матрица имеет размер только 3х3, вычисления еще более тривиальны.Однако я оставлю первую часть поста такой, как она есть, поскольку она описывает полностью правильное решение для общей матрицы любого размера, имеющего ранг 3.
Если ваша цель состоит в том, чтобы убить последнее единственное числоПри значении матрицы 3x3 svd на 3x3 кажется достаточно эффективным.Также будет некоторая потеря точности при генерации последнего сингулярного значения косвенными средствами.Например, сравните здесь наименьшее единственное значение, вычисленное с помощью svd, а затем с использованием собственных значений.Таким образом, вы можете увидеть здесь небольшую ошибку, поскольку формирование A '* A приведет к некоторой точности.Степень этой потери, вероятно, будет зависеть от номера условия A.
>> A = rand(3);
>> sqrt(eig(A'*A))
ans =
0.0138948003095069
0.080275195586312
1.50987693453097
>> svd(A)
ans =
1.50987693453097
0.0802751955863124
0.0138948003095054
Однако, если вы действительно хотите выполнить работу самостоятельно, вы можете это сделать.
- Вычислить характеристический полином непосредственно из матрицы 3x3 A '* A.
- Вычислить корни этого кубического полинома. Выберите самый маленький корень.
- Создайте соответствующий собственный вектор.
- Выполните обновление ранга один на А, убив часть А, которая лежит в направлении этого собственного вектора.
Является ли это проще или менее вычислительно эффективным, чем простой вызов svd с последующим обновлением ранга один? Я совсем не уверен, что это стоит усилий на 3х3. 3x3 SVD действительно очень быстро вычислить.