Die Nachahmung der Lösung der kleinsten Quadrate mit geringer Dichte in MATLAB führt zu möglichst genauen Ergebnissen in ⇐ Python

[Anonymous]

Ich habe einen Teil des MATLAB-Codes in Python umgeschrieben und nach viel Debugging bin ich an einem Punkt angelangt, an dem der Unterschied in den Ergebnissen zwischen MATLAB und Python auf die Ergebnisse der Lösung der kleinsten Quadrate zurückzuführen ist.
Die Ergebnisse scheinen sehr ähnlich.


Die Norm der Differenz (200-Elemente-Vektor) beträgt etwa 77,6. Die Norm der Residuen ist bis zur 7. Dezimalstelle gleich.
Aber leider unterscheidet sich das Ergebnis für mich erheblich. Das Endergebnis, das mein Code ausgibt, unterscheidet sich drastisch vom MATLAB-Ergebnis ... es sei denn, ich verwende MATLAB buchstäblich, um die Probleme der kleinsten Quadrate zu lösen. Ich habe mehrere verschiedene Löser für die kleinsten Quadrate von scipy.sparse.linalg ohne Erfolg ausprobiert.
Die Ergebnisse sind durchweg schlecht (und seltsamerweise sind die Gesamtergebnisse des gesamten Algorithmus auf die gleiche Weise durchweg falsch). Es ist fast so, als ob der Algorithmus ausschließlich für die Arbeit mit MATLAB erstellt wurde.
Im folgenden Code sind Lx,Ly Laplaceoperatoren eines einfachen Diagramms (daher PSD), Omega ist eine binäre Matrix.

Code: Select all

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