Имитирование разреженного решения методом наименьших квадратов MATLAB как можно точнее в Python ⇐ Python

[Anonymous]

Я переписал некоторый код MATLAB на Python и после долгой отладки достиг точки, когда разница в результатах между MATLAB и Python сводится к результатам решения методом наименьших квадратов.
Результаты кажется очень похожими.


Норма разницы (вектор из 200 элементов) составляет около 77,6. Норма остатков одинакова до 7-го знака после запятой.
Но, к сожалению для меня, результат существенно отличается. Конечный результат, который выдает мой код, радикально отличается от результата MATLAB... если только я буквально не использую MATLAB для решения задач наименьших квадратов. Я пробовал несколько разных решателей метода наименьших квадратов из scipy.sparse.linalg, но безрезультатно.
Результаты постоянно плохие (и, как ни странно, общие результаты всего алгоритма одинаково ошибочны). Это почти как если бы алгоритм был создан исключительно для работы с MATLAB.
В следующем коде Lx,Ly — лапласианы простого графа (следовательно, PSD), Omega — двоичная матрица.
from scipy import sparse
import numpy as np
import matlab

def solve_v(fixed_u, Lx, Ly, LxSquared, Omega, X_cols,RT_ju,_lambda,eng=None):
ninj = X_cols.shape[1]
uLy = Ly @ fixed_u
# (vLx @ vLx.T +2 (vLx @v) Ly ) + LySquared
y= RT_ju
Reg = LxSquared + (fixed_u.T @ uLy).flatten() * Lx + (uLy.T @ uLy).flatten()/(_lambda*4) * sparse.eye_array(Lx.shape[0])
B_L_0_5 = X_cols * np.sqrt(fixed_u.T @ (fixed_u * Omega))
if (eng is None):
x = sparse.linalg.cg(Reg,y)[0][:,np.newaxis]
else:
Reg_mat = matlab.double(Reg.toarray().tolist())
y_mat = matlab.double(y.tolist())
eng.workspace["Reg_mat"] = Reg_mat
eng.workspace["y_mat"] = y_mat
eng.eval("x = lsqminnorm(sparse(Reg_mat),y_mat);",nargout=0)
x = np.array(eng.workspace["x"])
if (eng is None):
BReg_inv_BL05 = np.empty((Ly.shape[0],ninj))
for i in range(B_L_0_5.shape[1]):
tmp =sparse.linalg.cg(Reg, B_L_0_5[:,i])
BReg_inv_BL05[:,i] = tmp[0]
else:
BL05_mat = matlab.double(B_L_0_5.tolist())
eng.workspace["BL05_mat"] = BL05_mat
eng.eval("BR_inv_B = lsqminnorm(sparse(Reg_mat),BL05_mat);",nargout=0)
BReg_inv_BL05 = np.array(eng.workspace["BR_inv_B"])
K = np.eye(ninj) + B_L_0_5.T @ BReg_inv_BL05
if (eng is None):
x = x-BReg_inv_BL05 @ np.linalg.lstsq(K,BReg_inv_BL05.T @ y)[0]
else:
K_mat = matlab.double(K)
eng.workspace["K"] = K_mat
eng.workspace["y_mat"] = matlab.double(BReg_inv_BL05.T @ y)
eng.eval("tmp = lsqminnorm(K,y_mat);",nargout=0)
x = x - BReg_inv_BL05 @ np.array(eng.workspace["tmp"])

return x


Подробнее здесь: https://stackoverflow.com/questions/798 ... s-possible