NumPy 中最小二乘算法的高效计算
Efficient computation of the least-squares algorithm in NumPy
我需要在最小二乘意义上求解一大组线性系统。我无法理解 numpy.linalg.lstsq(a, b)
、np.dot(np.linalg.pinv(a), b)
和数学实现的计算效率差异。
我使用以下矩阵:
h=np.random.random((50000,100))
a=h[:,:-1].copy()
b=-h[:,-1].copy()
算法的结果是:
# mathematical implementation
%%timeit
np.dot(np.dot(np.linalg.inv(np.dot(a.T,a)),a.T),b)
10 次循环,3 次循环中的最佳:36.3 毫秒 每个循环
# numpy.linalg.lstsq implementation
%%timeit
np.linalg.lstsq(a, b)[0]
10 次循环,3 次中的最佳循环:103 毫秒 每个循环
%%timeit
np.dot(np.linalg.pinv(a), b)
1 个循环,3 个循环中的最佳:每个循环 216 毫秒
为什么会有差异?
例程 lstsq
处理任何系统:over-determined、under-determined 或 well-determined。它的输出是你从 pinv(a)*b 得到的,但它比计算伪逆更快。原因如下:
一般建议:除非您确实需要,否则不要计算逆矩阵。为特定的右手边求解系统比反转其矩阵更快。
然而,你通过求解 aTa = aTb 的方法更快,即使你正在反转矩阵。是什么赋予了?问题是,反转 aTa 仅在 a 具有完整列等级时才有效。因此,您已将问题限制在这种特定情况下,并提高了速度,作为对通用性较低以及安全性较低的权衡,如下所示。
但是求逆矩阵仍然效率低下。如果您知道 a 具有完整的列排名,则以下操作比您的三次尝试中的任何一次都快:
np.linalg.solve(np.dot(a.T, a), np.dot(a.T, b))
也就是说,在处理条件较差的矩阵时,lstsq
仍然比上面的更好。形成乘积 aTa 基本上是条件数的平方,所以你更有可能得到无意义的结果。这是一个警示性的例子,使用了 SciPy 的 linalg 模块(本质上等同于 NumPy,但有更多的方法):
import numpy as np
import scipy.linalg as sl
a = sl.hilbert(10) # a poorly conditioned square matrix of size 10
b = np.arange(10) # right hand side
sol1 = sl.solve(a, b)
sol2 = sl.lstsq(a, b)[0]
sol3 = sl.solve(np.dot(a.T, a), np.dot(a.T, b))
这里lstsq
给出了与solve
几乎相同的输出(这个系统的唯一解)。然而,sol3
由于数字问题而完全错误(你甚至不会被警告)。
sol1:
[ -9.89821788e+02, 9.70047434e+04, -2.30439738e+06,
2.30601241e+07, -1.19805858e+08, 3.55637424e+08,
-6.25523002e+08, 6.44058066e+08, -3.58346765e+08,
8.31333426e+07]
sol2:
[ -9.89864366e+02, 9.70082635e+04, -2.30446978e+06,
2.30607638e+07, -1.19808838e+08, 3.55645452e+08,
-6.25535946e+08, 6.44070387e+08, -3.58353147e+08,
8.31347297e+07]
sol3:
[ 1.06913852e+03, -4.61691763e+04, 4.83968833e+05,
-2.08929571e+06, 4.55280530e+06, -5.88433943e+06,
5.92025910e+06, -5.56507455e+06, 3.62262620e+06,
-9.94523917e+05]
我需要在最小二乘意义上求解一大组线性系统。我无法理解 numpy.linalg.lstsq(a, b)
、np.dot(np.linalg.pinv(a), b)
和数学实现的计算效率差异。
我使用以下矩阵:
h=np.random.random((50000,100))
a=h[:,:-1].copy()
b=-h[:,-1].copy()
算法的结果是:
# mathematical implementation
%%timeit
np.dot(np.dot(np.linalg.inv(np.dot(a.T,a)),a.T),b)
10 次循环,3 次循环中的最佳:36.3 毫秒 每个循环
# numpy.linalg.lstsq implementation
%%timeit
np.linalg.lstsq(a, b)[0]
10 次循环,3 次中的最佳循环:103 毫秒 每个循环
%%timeit
np.dot(np.linalg.pinv(a), b)
1 个循环,3 个循环中的最佳:每个循环 216 毫秒
为什么会有差异?
例程 lstsq
处理任何系统:over-determined、under-determined 或 well-determined。它的输出是你从 pinv(a)*b 得到的,但它比计算伪逆更快。原因如下:
一般建议:除非您确实需要,否则不要计算逆矩阵。为特定的右手边求解系统比反转其矩阵更快。
然而,你通过求解 aTa = aTb 的方法更快,即使你正在反转矩阵。是什么赋予了?问题是,反转 aTa 仅在 a 具有完整列等级时才有效。因此,您已将问题限制在这种特定情况下,并提高了速度,作为对通用性较低以及安全性较低的权衡,如下所示。
但是求逆矩阵仍然效率低下。如果您知道 a 具有完整的列排名,则以下操作比您的三次尝试中的任何一次都快:
np.linalg.solve(np.dot(a.T, a), np.dot(a.T, b))
也就是说,在处理条件较差的矩阵时,lstsq
仍然比上面的更好。形成乘积 aTa 基本上是条件数的平方,所以你更有可能得到无意义的结果。这是一个警示性的例子,使用了 SciPy 的 linalg 模块(本质上等同于 NumPy,但有更多的方法):
import numpy as np
import scipy.linalg as sl
a = sl.hilbert(10) # a poorly conditioned square matrix of size 10
b = np.arange(10) # right hand side
sol1 = sl.solve(a, b)
sol2 = sl.lstsq(a, b)[0]
sol3 = sl.solve(np.dot(a.T, a), np.dot(a.T, b))
这里lstsq
给出了与solve
几乎相同的输出(这个系统的唯一解)。然而,sol3
由于数字问题而完全错误(你甚至不会被警告)。
sol1:
[ -9.89821788e+02, 9.70047434e+04, -2.30439738e+06,
2.30601241e+07, -1.19805858e+08, 3.55637424e+08,
-6.25523002e+08, 6.44058066e+08, -3.58346765e+08,
8.31333426e+07]
sol2:
[ -9.89864366e+02, 9.70082635e+04, -2.30446978e+06,
2.30607638e+07, -1.19808838e+08, 3.55645452e+08,
-6.25535946e+08, 6.44070387e+08, -3.58353147e+08,
8.31347297e+07]
sol3:
[ 1.06913852e+03, -4.61691763e+04, 4.83968833e+05,
-2.08929571e+06, 4.55280530e+06, -5.88433943e+06,
5.92025910e+06, -5.56507455e+06, 3.62262620e+06,
-9.94523917e+05]