时间戳

首页 > 解决方案 > Numpy/Scipy：用相同的设计矩阵求解几个最小二乘

numpy - Numpy/Scipy：用相同的设计矩阵求解几个最小二乘

问题描述

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

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

想法？

标签： numpyscipyleast-squares

解决方案

如果您的线性系统已确定，我将存储MLU 分解并将其单独用于所有' 或仅对表示水平堆叠' 的b2d 数组进行一次求解调用，这实际上取决于您的问题但这是全球相同的想法。假设您一次获得了每一个，那么：Bbb

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)])

M_pinv或者简单地使用pinv（基于 lstsq）或pinv2（基于 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例如，只需存储 的 SVD M，然后手动解决：

# 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

不过，您可以使用适当的尺寸检查这些内容。

希望这可以帮助。

推荐阅读