python - 如何为向量输入和观测值建立线性最小二乘回归模型

Question

假设我们有一个由矩阵相关的输入和观察系统:

如果我们有一组观测值 y，基于一组输入 x，我可以设置一个非线性最小二乘例程来拟合矩阵 m 的参数

%pylab inline
from scipy.optimize import least_squares
n_observations = 100
m = random.random(16).reshape(4, 4)
x = random.random(n_observations*4).reshape(n_observations, 4, 1)
noise = (random.random(n_observations*4).reshape(n_observations, 4, 1)-0.5) * 0.01
y = einsum('ij,njk->nik', m, x)

def residuals(x0):
    return (y + noise - einsum('ij,njk->nik', x0.reshape(4, 4), x)).flatten()

res = least_squares(residuals, x0=random.random(16))
m_fit = res.x.reshape(4, 4)
diff = m_fit - m
print('     m actual | m fit    | diff     ')
print('     -------- | -------- | ---------')
for i in range(4):
    for j in range(4):
        print(f'm{i+1}{j+1}: {m[i,j]:0.06f} | {m_fit[i,j]:0.06f} | {diff[i,j]:+0.06f}')

>>> (for example)
     m actual | m fit    | diff
     -------- | -------- | --------
m11: 0.259722 | 0.259461 | -0.000261
m12: 0.266986 | 0.266999 | +0.000012
m13: 0.373180 | 0.373662 | +0.000482
m14: 0.570387 | 0.569813 | -0.000574
m21: 0.462023 | 0.462099 | +0.000076
m22: 0.875758 | 0.876651 | +0.000893
m23: 0.420369 | 0.419884 | -0.000485
m24: 0.335546 | 0.334505 | -0.001041
m31: 0.625779 | 0.626269 | +0.000490
m32: 0.499375 | 0.499400 | +0.000025
m33: 0.871075 | 0.870183 | -0.000892
m34: 0.497999 | 0.498878 | +0.000879
m41: 0.367814 | 0.366537 | -0.001277
m42: 0.020419 | 0.020412 | -0.000007
m43: 0.221916 | 0.221764 | -0.000153
m44: 0.758361 | 0.759409 | +0.001048

我的问题是，是否可以使用线性最小二乘回归方法(a.la)来做到这一点。 numpy.linalg.lstsq

我对线性回归不是很熟悉，但看起来这应该是可能的，我只是不知道如何设置问题来以这种方式执行它。 numpy.linalg.lstsq 方法似乎不是为处理这种特定场景而设置的，我不确定还有什么可以，所以我正在寻找一些关于这方面的指导。

Answer 1

编写方程的另一种方法是 X*M = Y，其中 X 是 (n_obs, 4) 输入矩阵，M 是(4, 4) 未知数矩阵，Y 是 (n_obs, 4) 观测值。

然后numpy.linalg.lstsq可以使用:

import numpy as np
from scipy.optimize import least_squares

n_observations = 100

np.random.seed(seed=1234)
X = np.random.random((n_observations, 4))
M = np.random.random((4, 4))

Y = np.einsum('ni,ij->nj', X, M)

M_fit, residuals, rank, s = np.linalg.lstsq(X, Y, rcond=None)

print(M)
#[[0.71499388 0.72409148 0.01867644 0.2858131 ]
# [0.58048634 0.93078663 0.3389969  0.12008312]
# [0.51627271 0.69920706 0.29864068 0.86160962]
# [0.9058072  0.76858325 0.26123164 0.9384556 ]]

print(M_fit)
#[[0.71499388 0.72409148 0.01867644 0.2858131 ]
# [0.58048634 0.93078663 0.3389969  0.12008312]
# [0.51627271 0.69920706 0.29864068 0.86160962]
# [0.9058072  0.76858325 0.26123164 0.9384556 ]]

通过编写方程的转置，即 M' * X' = Y'，可以检索您的符号。