Итак, я пытаюсь найти обобщенную обратную A+ матрицы A, используя svd
это мой код для SVD
A4 = np.array([[3, 1,2,3,4],
[5,2,3,4,5],
[7,4,5,6,7]])
print(A4)
U4, s4, V4 = np.linalg.svd(A4)
D4 = sp.linalg.diagsvd(s4, 3, 5)
print('UDVt =\n', U4@D4@V4 )
Который отлично работает и дает мне правильную матрицу
И это код, который я пытаюсь использовать, чтобы получить обобщенное обратное
D4p = np.array([[D4[0,0]**-1,0,0],
[0,D4[1,1]**-1,0],
[0,0,D4[2,2]**-1],
[0,0,0],
[0,0,0]])
A4plus = np.transpose(V4) @ D4p @ np.transpose(U4)
print(np.round(A4plus @ A4,2))
print(np.round(A4 @ A4plus,2))
Я попытался изменить транспонированную матрицу.
Разве не этого вы ожидали?
@chrslg Да, это то, чего я ожидал от A4@A4plus, но A4plus@A4 неправильный. Но я согласен с вами, что он дает то же самое, что и np.linalg.pinv(A4). SVD для меня новичок, поэтому мне было интересно, почему результаты такие разные, и я подумал, что допустил ошибку в коде. Но, кажется, с кодом все в порядке.
Не ждите, что A4plus @ A4
будет идентификационной матрицей. Обычно этого достичь невозможно, так как у вас недостаточно ранга строки. А A4 @ A4plus
— это идентичность, как и должно быть.
Предполагая, что строки не являются линейно зависимыми (иначе S-матрица имеет элементы нулевого значения и я не смогу их инвертировать), можно сказать
U * S * Vt = A
=> A * (Vt.T * 1/S * U.T) = I
В коде:
from numpy import array, round
from numpy.linalg import svd
from scipy.linalg import diagsvd
A = array([[3, 1, 2, 3, 4],
[5, 2, 3, 4, 5],
[7, 4, 5, 6, 7]])
U, S, Vt = svd(A)
S_inv = diagsvd(1/S, A.shape[1], A.shape[0])
A_inv = Vt.T @ S_inv @ U.T
print(round(A @ A_inv, 5))
В результате это дает псевдообратное A
[[ 1. 0. -0.] [ 0. 1. 0.] [-0. -0. 1.]]
Вот соответствующий абзац по этому поводу: https://inst.eecs.berkeley.edu/~ee127/sp21/livebook/def_pseudo_inv.html
Хорошо, последний комментарий проясняет, в чем заключается настоящий вопрос: не о SVD как таковом, а о псевдообратном. Почему, если A⁺ является псевдоинверсией, A⁺A не является идентифицирующей (хотя AA⁺ является)?
Не имеет значения, как вычисляется A⁺.
Ap = np.linalg.pinv(A)
Ap = [email protected]([email protected])
U,s,V = np.linalg.svd(A); Ap = ([email protected](1/s, 3, 5)@V).T
так как результат тот же.
Итак, ответ на реальный вопрос прост: псевдоинверсия не является инверсией. Мы используем его, потому что обратного не существует.
Если бы A⁺A также было тождественным, то A⁺ было бы инверсией A. ∀u,v, Au=v ⇔ A⁺v=u, будет означать, что A обратимо, а его обратным является A⁺.
Чего не может быть. Так как это матрица 3×5.
Думайте об этом как о теории информации. Если бы A⁺A было тождеством, это означало бы, что если у вас есть 5 чисел для хранения в файле, вы могли бы поместить их в вектор U=[u₁,u₂,u₃,u₄,u₅]. Вычислите V=AV. Это дает вам только 3 числа V=[v₁,v₂,v₃]. Сохраните эти три числа в файле. Затем, когда вы захотите получить файл, просто вычислите U'=A⁺V, что в точности равно U.
Идеальный алгоритм сжатия. Это будет работать для всех данных U. ∀U, V=AU, тогда U'=A⁺V ⇒ U'=A⁺AU = IU = U. Не может быть правдой. Вы не можете представить весь возможный набор из 5 чисел, из которых 3 числа. Вы не можете вместить все 5D-пространство в 3D-пространство.
Алгоритм сжатия, конечно, существует. Но они работают только для подмножества возможных данных. Они работают, потому что мы знаем, что фактические данные, которые люди на самом деле используют, представляют собой лишь очень-очень небольшое подмножество всех возможных данных. Сжатие изображений работает только потому, что 99,99999999% изображений выглядят для нас как белый шум, и только оставшиеся 0,00000001% выглядят как изображение чего-либо. В противном случае это было бы невозможно.
Возвращаясь в линейный мир, вот что у вас здесь есть. «Сжатие» 5 чисел до 3, которое работает только для 3D-подмножества вашего 5D-пространства. То есть, если ваши данные представляют собой линейную комбинацию строк A.
Например, для U=[1,2,3,4,5] AU = [51,69,103]. И A⁺AU = [1,2,3,4,5]. Но это потому, что [1,2,3,4,5] — это 4*A₁-5A₂+2A₃ (3 строки A).
Итак, если вы случайно знаете, что U, с которым вы сталкиваетесь, лежат не во всем пятимерном пространстве, а только в небольшом его трехмерном подпространстве, то вы действительно можете ожидать, что эта псевдообратная работа будет работать (то есть, кстати, СВД - это нахождение действительно актуального подпространства).
Но если U может быть чем угодно, а не просто линейной комбинацией строк A, например [1,0,4,8,2]. И AU — это [43,59,89], и тогда A⁺AU — это [1,2,3,4,5], а не [1,0,4,8,2] (конечно, я выбрал свой Например, A@[1,2,3,4,5] = A@[1,0,4,8,2]. Но дело в том, что оно не может быть уникальным. Вы не можете иметь биекцию между 3D-пространствами. и 5D-пространство).
Кстати, [1,2,3,4,5] = A⁺@A@[1,0,4,8,2] не является каким-либо другим u, например Au=A[1,0,4,8,2 ]; а не просто любой другой вектор, который разделяет результат Au с [1,0,4,8,2]. Это тот, который представляет собой линейную комбинацию строк А. А это означает другое. Это означает, что A⁺A[1,0,4,8,2] — это линейная комбинация строк A, ближайшая (с точки зрения квадратической разности) к [1,0,4,8,2].
Тот факт, что A⁺A@[1,0,4,8,2] = [1,2,3,4,5] означает, что [1,2,3,4,5] является линейной комбинацией строк A, которые минимизируют ||AU-[1,0,4,8,2]||
Это еще один способ увидеть псевдоинверсию.
Если бы A была матрицей ранга 5. Тогда линейные комбинации строк A будут составлять все 5D-пространство. Другими словами, любой вектор X будет линейной комбинацией строк A. Итак, вы ожидаете, что A⁺AX будет X. Потому что X является линейной комбинацией строк A, и среди тех, к которым можно быть ближе X, чем сам X?
Итак, A⁺AX=X ∀X. Итак, A⁺A=Id. Таким образом, A⁺ — это не просто псевдоинверсия. Это обратное.
Но поскольку A не имеет ранга 5, то не все векторы X представляют собой линейную комбинацию строк A. И затем, если вы хотите найти линейную комбинацию строк A, напоминающих X, вам нужно решить задачу наименьших квадратов и найти Y, например, Y представляет собой линейную комбинацию строк A и ||X-Y|| является минимальным. Оно не может быть 0 каждый раз. Но вы можете подойти как можно ближе. Сказав Y=A⁺AX
Или, используя комбинации строк: A⁺ᵀ X — это коэффициенты строк A, наиболее близких к X. Например, для X=[1,0,4,8,2] Ap.T@X = [4 ,-5,2], что означает, что вы не можете получить αA₁+βA₂+γA₃ ближе к X, чем при α=4, β=-5 γ=2. И тогда αA₁+βA₂+γA₃ = A⁺AX =[1,2,3,4,5].
Но, конечно, если бы это было [1,0,4,8,2] вместо [1,2,3,4,5], это означало бы, что любой вектор X представляет собой линейную комбинацию строк A. И это не может быть, да? Никакой 5D-вектор не может быть линейной комбинацией трех векторов.
A⁺A не является тождеством, потому что в противном случае оно было бы обратным, а не псевдоинверсным.
AA⁺ является тождественным, потому что AA⁺Y (Y представляет собой 3D-вектор) является ближайшей к Y линейной комбинацией столбцов A. Если A имеет ранг 3, существует одна линейная комбинация столбцов A, которая в точности соответствует Y для всех Y. Следовательно, AA⁺=Id
A⁺A не является тождественным, потому что A⁺AX (X представляет собой 5D-вектор) является ближайшей к X линейной комбинацией строк A. И вы не можете иметь точную комбинацию трех строк A, которая является X, для всех X.
Так? Это странный способ сделать это (со всем этим ручным расширением). Но, кажется, это работает, не так ли? Вы получаете именно то, что дает
np.linalg.pinv(A4)
(я не знаю, делает ли он то же самое или нет. SVD — это один из способов вычисления обратного Мура-Пенроуза, но не единственный). Или именно то, чтоA4.T @ np.linalg.inv([email protected])
дает. Это другой путь. А A4@A4plus — это более или менее идентичность.