Что происходит

Вы пытаетесь вычислить псевдообратную Мура-Пенроуза матрицы, ранг которой равен 1, а размерность намного больше 1.

Это законно, конечно. Это даже то, для чего псевдоинверсия.

Но в какой-то момент вычислений будет задействовано множество малых значений, которые должны быть точно равны 0, но в числовом выражении это не так.

Давайте сделаем это вручную

import numpy as np

AA=np.asarray([[0.1,0.1]]*100)

u,d,vh=np.linalg.svd(AA, full_matrices=False)

# d are the singular value
# And AA is (u*d)@vh
print((u*d)@vh)

Дает

# [[0.1, 0.1],
#  [0.1, 0.1],
#  ... (97 lines)
#  [0.1, 0.1]]

Тогда псевдоинверсия

InvAA = (vh.T*(1/d)) @ u.T

Ну, почти так. Потому что d содержит некоторый «псевдо 0». Поэтому мы должны использовать псевдоинверсию d, чтобы выполнить это вычисление. Псевдообратная диагональная матрица проста: это матрица, диагональ которой состоит из обратной ненулевых диагональных элементов.

Здесь d — это 1D и просто значения диагонали. Так что invd будет [1/x for x in d if np.abs(x)>0]

Но, поскольку это численное вычисление, 0 — идеальный порог. Нам нужен более реалистичный. Скажем 1e-15

В нашем случае с этим порогом все в порядке: мы знаем, что истинный ранг равен 1. И вычисление svd дает для d: d=[1.41421356e+00, 1.79973556e-16]

(Таким образом, это действовало так, как если бы ранг был равен 2 с очень маленьким сингулярным значением)

invd=[1/x for x in d if np.abs(x)>1e-15]
InvAA = (vh.T * invd) @ u.T
print(InvAA)

Вы можете проверить, что это работает:
InvAA@AA@InvAA это InvAA
AA@InvAA@AA это AA
AA@InvAA и InvAA@AA оба симметричны

Давайте сделаем то же самое из AA.T

u,d,vh=np.linalg.svd(AA.T, full_matrices=False)
invd=[1/x for x in d if np.abs(x)>1e-15]
InvAAT = (vh.T * invd) @ u.T
print(AA.T@[email protected])

На этот раз это не работает. И если мы посмотрим на значения InvAAT, то произойдет очевидное: очень большие значения. Обратите внимание, что это почти работает. AA.T@[email protected] не так уж далеко от AA.T по сравнению с огромными значениями, из которых состоят вычисления. Таким образом, это, очевидно, проблема числовой точности.

Если мы посмотрим на d, мы сразу поймем, почему

d это [1.41421356e+00, 1.95277784e-15]

Итак, нашего порога 1e-16 здесь недостаточно.

Что, если мы просто понизим его

u,d,vh=np.linalg.svd(AA.T, full_matrices=False)
invd=[1/x for x in d if np.abs(x)>1e-14]
InvAAT = (vh.T * invd) @ u.T
print(AA.T@[email protected])

Он снова работает!

Конечно, с риском потери информации, если одно из этих вырезанных «почти 0» сингулярных значений не было бы на самом деле 0. Но в нашем случае нет проблем (в любом случае, мы знаем, что реальный ранг равен 1. Так что мы могли бы в действительности выберите просто сохранить первое значение d и отбросить остальные).

Почему так происходит

Теперь, почему это чаще происходит с AA.T, чем с AA; почему числовая ошибка больше с AA.T, чем с AA. Я не могу сказать, так как я не знаю специфики алгоритма svd numpy.

Но если вы посмотрите на множество эквивалентных способов вычисления псевдообратного метода Мура-Пенроуза или разложения по сингулярным значениям, вы увидите, что он включает в себя либо вычисление собственных значений M.T@M. А AA.T@AA — это матрица 2×2. Когда [email protected] - матрица размера nn × nn. Так что неудивительно, что последний включает больше шума.

Другой метод - вычислить предел (M.T@M+εI)⁻¹@M.T. Что также включает инвертирование матрицы 2 × 2 с рангом почти 1, но с одним значением ε, чтобы сделать ранг равным 2 в случае AA. И матрица 100 × 100 (для nn = 100, как в моем примере) с рангом почти 1, но 99 значений ε для ранга 100 в случае AA.T. Очевидно, вам будет труднее заставить его сходиться, и числовая ошибка, когда ε действительно мала, в последнем случае будет больше, чем в первом.

Вы можете играть с этим методом

np.linalg.inv(AA.T@AA + 1e-13*np.identity(2))@AA.T

Кажется, дает разумное значение для псевдообратного. Пока

np.linalg.inv(AA.T@AA + 1e-14*np.identity(2))@AA.T

хуже (поэтому, если вы делаете ε все меньше и меньше, аппроксимация становится все лучше и лучше до 1e-13, затем она становится все хуже и хуже. По обычной причине: математика говорит «чем меньше, тем лучше», а числовая ошибка говорит «чем меньше, тем лучше»). чем меньше, тем хуже». Судя по всему, числовая ошибка лидирует в районе 1e-14.

Если вы сделаете то же самое для случая AA.T

np.linalg.inv(AA.T@AA + ε*np.identity(2))@AA.T

На этот раз это примерно 1e-10, когда числовая ошибка начинает ухудшать результаты.

Итак, с этим также вы можете видеть, что с численной точки зрения сложнее вычислить псевдоинверсию AA.T, чем вычислить псевдоинверсию AA

Ни один из двух приведенных мной методов, вероятно, не используется numpy. Но такое же соображение, вероятно, применимо. Числовая ошибка просто больше в случае AA.T

Как это предотвратить

Итак, короче говоря: то, что вы должны считать 0 при выборе сингулярных значений, должно быть больше во втором случае. А с pinv именно для этого и нужен параметр rcond.

Просто добавьте rcond=1e-10 к звонку pinv, и вы получите это

Но, конечно, это связано с обычным компромиссом: слишком большой rcond риск игнорирует значительные (маленькие, но на самом деле не 0, а не просто числовую ошибку) значения. Таким образом, вы рискуете получить псевдоинверсию ранга, даже меньшего, чем ранг вашей матрицы. Конечно, не в вашем случае, так как ваша матрица имеет ранг 1. И, кроме нулевой матрицы, меньшего ранга в любом случае нет. Как я уже сказал, вы можете даже просто с помощью метода svd удалить все значения d, кроме первого: (v.h/d[0])@u. Но это только потому, что вы знаете, что ваша матрица действительно имеет ранг 1 (и сингулярное значение сортировки svd numpy от большего к меньшему)