Numpy/Scipy : 用相同的设计矩阵求解几个最小二乘法

Question

我遇到了一个通过 scipy.linalg.lstsq(M,b) 解决的最小二乘问题，其中：

M 的形状为 (n,n)
b 的形状为 (n,)

问题是我必须花很多时间来解决不同的 b。我怎样才能做更有效率的事情？我猜 lstsq 做了很多独立于 b 的值的事情。

想法？

Answer 1

您的问题不清楚，但我猜您的意思是针对不同的数组 (b0, b1, b2..) 计算方程式 Mx=b 到 scipy.linalg.lstsq(M,b)。如果是这种情况，您可以将过程与 concurrent.futures.ProcessPoolExecutor 并行化。这方面的文档相当简单，可以同时帮助 python 运行多个 scipy 求解器。

希望对您有所帮助。

Answer 2

您可以将 M 分解为 QR 或 SVD 产品并手动找到 lsq 解决方案。

Answer 3

如果您的线性系统已确定，我将存储 M LU 分解并将其用于所有 b 单独或简单地为 2d- 执行一次求解调用数组 B 表示水平堆叠的 b，这实际上取决于您的问题，但这在全球范围内都是相同的想法。假设您一次得到每个 b 个，那么：

import numpy as np
from scipy.linalg import lstsq, lu_factor, lu_solve, svd, pinv

# as you didn't specified any practical dimensions
n = 100
# number of b's
nb_b = 10

# generate random n-square matrix M
M = np.random.rand(n**2).reshape(n,n)

# Set of nb_b of right hand side vector b as columns
B = np.random.rand(n*nb_b).reshape(n,nb_b)

# compute pivoted LU decomposition of M
M_LU = lu_factor(M)
# then solve for each b
X_LU = np.asarray([lu_solve(M_LU,B[:,i]) for i in range(nb_b)])

但如果是欠定或多定，你需要像你一样使用lstsq：

X_lstsq = np.asarray([lstsq(M,B[:,i])[0] for i in range(nb_b)])

或者简单地用 pinv (built on lstsq) or pinv2 存储伪逆 M_pinv（基于 SVD）：

# compute the pseudo-inverse of M
M_pinv = pinv(M)
X_pinv = np.asarray([np.dot(M_pinv,B[:,i]) for i in range(nb_b)])

或者你也可以自己做，比如pinv2，只存储M的SVD，然后手动解决这个问题：

# compute svd of M
U,s,Vh = svd(M)

def solve_svd(U,s,Vh,b):
    # U diag(s) Vh x = b <=> diag(s) Vh x = U.T b = c
    c = np.dot(U.T,b)
    # diag(s) Vh x = c <=> Vh x = diag(1/s) c = w (trivial inversion of a diagonal matrix)
    w = np.dot(np.diag(1/s),c)
    # Vh x = w <=> x = Vh.H w (where .H stands for hermitian = conjugate transpose)
    x = np.dot(Vh.conj().T,w)
    return x

X_svd = np.asarray([solve_svd(U,s,Vh,B[:,i]) for i in range(nb_b)])

如果使用 np.allclose 检查，它们都会给出相同的结果（除非系统没有很好地确定导致 LU 直接接近失败）。最后在性能方面：

%timeit M_LU = lu_factor(M); X_LU = np.asarray([lu_solve(M_LU,B[:,i]) for i in range(nb_b)])
1000 loops, best of 3: 1.01 ms per loop

%timeit X_lstsq = np.asarray([lstsq(M,B[:,i])[0] for i in range(nb_b)])
10 loops, best of 3: 47.8 ms per loop

%timeit M_pinv = pinv(M); X_pinv = np.asarray([np.dot(M_pinv,B[:,i]) for i in range(nb_b)])
100 loops, best of 3: 8.64 ms per loop

%timeit U,s,Vh = svd(M); X_svd = np.asarray([solve_svd(U,s,Vh,B[:,i]) for i in range(nb_b)])
100 loops, best of 3: 5.68 ms per loop

尽管如此，您还是要用适当的尺寸来检查这些。

希望这对您有所帮助。

Numpy/Scipy : 用相同的设计矩阵求解几个最小二乘法

Numpy/Scipy : solving several least squares with the same design matrix

numpy

scipy

least-squares