求解 numpy/scipy 中的 3D 最小二乘法

Question

对于一些 100 左右的整数 K，我有 2 * K(n, n)数组：X_1, ..., X_K和Y_1, ..., Y_K。

我想同时执行 K 最小二乘法，即求 n 乘 n 矩阵A最小化 k：sum_k norm(Y_k - A.dot(X_k), ord='fro') ** 2上的平方和(A不得依赖于k(。

我正在寻找一种简单的方法来使用 numpy 或 scipy 来做到这一点。 我知道我想要最小化的函数是 A 中的二次形式，所以我可以手动完成，但我正在寻找一种现成的方法。有吗？

Answer 1

如果n是一个小数字，这样的东西就可以工作。

import numpy as np
from scipy.optimize import minimize
K = 5
n = 10
X = np.random.random_sample((K, n, n))
Y = np.random.random_sample((K, n, n))
def opt(A):
A = np.reshape(A, (n, n))
# maybe need to transpose X.dot(a) ?
# if axis is a 2-tuple, it specifies the axes that hold 2-D matrices, 
# and the matrix norms of these matrices are computed.
return np.sum(np.linalg.norm(Y - X.dot(A), ord='fro', axis=(1, 2)) ** 2.0)
A_init = np.random.random_sample((n, n))
print(minimize(opt, A_init ))

注意：minimize默认使用的优化算法是本地的。

Answer 2

事实上答案很简单，我只需要通过水平堆叠Y_k(创建 Y(和X_k(创建 X(来创建更大的矩阵 Y 和 X。然后我可以解决一个常规的 2d 最小二乘问题：最小化norm(Y - A.dot(X))

Answer 3

我不能帮助python，但这是数学解决方案，以防它有帮助。 我们力求最小化

E = Sum { Tr (Y[j]-A*X[j])*(Y[j]-A*X[j])'}

一些代数产量

E = Tr(P-A*Q'-Q*A'+A*R*A')
where
P = Sum{ Y[j]*Y[j]'}
Q = Sum{ Y[j]*X[j]'}
R = Sum{ X[j]*X[j]'}

如果 R 是可逆的，则多一点代数产生

E = Tr( (A-Q*S)*R*(A-Q*S)') + Tr( P - Q*S*Q')
where S = inv( R)

因为

(A-Q*S)*R*(A-Q*S)' is positive definite, 

我们通过取 A = Q*S 来最小化 E。

在这种情况下，算法将是：

compute Q
compute R
solve A*R = Q for A (eg by finding the cholesky factors of R)

如果 R 不可逆，我们应该对 S 使用广义逆而不是普通逆。