SVD-разложение для решения линейных уравнений ⇐ Python

[Anonymous]

Я ищу здесь Algo и пытаюсь реализовать этот код, получаю разные l2-нормы для результирующих векторов параметров для линейного уравнения.
Где я ошибаюсь, пытаясь внедрить код?

Код: Выделить всё

import numpy as np
from scipy import linalg

np.random.seed(123)
v = np.random.rand(4)
A = v[:,None] * v[None,:]
b = np.random.randn(4)

x = linalg.inv(A.T.dot(A)).dot(A.T).dot(b) #Usually not recommended because of Numerical Instability of the Normal Equations  https://johnwlambert.github.io/least-squares/
l2_0= linalg.norm(A.dot(x) - b)
print("manually: ", l2_0)

x = linalg.lstsq(A, b)[0]
l2_1= linalg.norm(A.dot(x) - b)
print("scipy.linalg.lstsq: ", l2_1)

# 2-norm of two calculations compared
print(np.allclose(l2_0, l2_1, rtol=1.3e-1))

def direct_ls_svd(x,y):
# append a columns of 1s (these are the biases)
x = np.column_stack([np.ones(x.shape[0]), x])

# calculate the economy SVD for the data matrix x
U,S,Vt = linalg.svd(x, full_matrices=False)

# solve Ax = b for the best possible approximate solution in terms of least squares
x_hat = Vt.T @ linalg.inv(np.diag(S)) @ U.T @ y
#print(x_hat)

# perform train and test inference
#y_pred = x @ x_hat

return y-x @ x_hat     #x_hat

x= direct_ls_svd(A, b)
l2_svd= linalg.norm(A.dot(x) - b)
print("svd: ", l2_svd)

# LU
x= linalg.solve(A.T@A, A.T@b)
l2_solve= linalg.norm(A.dot(x) - b)
print("scipy.linalg.solve: ", l2_solve)

# manually:  2.9751344995811313
# scipy.linalg.lstsq:  2.9286130558050654
# True
# svd:  6.830550019041984
# scipy.linalg.solve:  2.928613055805065

Если моя ошибка связана с реализацией алгоритма декомпозиции SVD для решения задачи наименьших квадратов или, возможно, с относительным округлением Numpy Scipy или различиями в точности? Как исправить svd-алгоритм для наименьших квадратов, чтобы он стал сопоставимым с алгоритмом scipy? И будет ли этот алгоритм быстрее и меньше потреблять памяти, чем итеративные методы наименьших квадратов?
PS:
Приложения SVD или здесь, SVD для PCA и PLS-SVD - моя конечная цель - будет ли алгоритм таким же, как и для аппроксимации наименьших квадратов? Меня вообще смущает такой вопрос (даже с примерами кода). Может ли кто-нибудь добавить немного ясности для новичков вроде меня?
PPPS: применяя такую ​​реализацию - я получаю еще худший результат: svd: 10.031259300735462 для l2-norm
PPPS: также мне не хватает некоторого понимания svd в сингулярном спектральном разложении, если существует ссылка, как 1-й для неконтролируемого dim.reduction и 2-й для непараметрического анализа TS, для практики
PPPPS: ! PCA используется для предварительной оценки, если существует мультиколлинеарность, в противном случае могут быть получены странные результаты (предвзятость и т. д.)... (если нет коллинеарности => нет чувствительности к ошибке оценки, то есть показано в небольшом количестве условий анализа OLS, - vc.vs ​​огромный cond.num.== коллинеарность для многомерной/многомерной регрессии)

Подробнее здесь: https://stackoverflow.com/questions/776 ... on-solving