在 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

我对 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 的广义逆而不是普通逆。

Answer 2

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

Answer 3

如果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 默认使用的优化算法是本地的。