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

Question

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

我想同时执行 K 最小二乘法，即找到 n x 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

我无法帮助 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 3

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

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

Solve 3D least squares in numpy/scipy

python

numpy

linear-algebra

scipy

least-squares