# NumpyとScipyでlinalg.lstsqを比較する

numpy tutorialの一環としてnumpy.linalg.lstsqをやる。この関数についてはこのサイトに以下のように書いてある。ついでに、scipy版とも何が違うのか比較してみる。

Solves the equation $ax = b$ by computing a vector $x$ that minimizes the Euclidean 2-norm $|| b – ax ||^2$. The equation may be under-, well-, or over- determined (i.e., the number of linearly independent rows of a can be less than, equal to, or greater than its number of linearly independent columns). If a is square and of full rank, then $x$ (but for round-off error) is the “exact” solution of the equation.

スポンサーリンク

## numpy.linalg.lstsq¶

Fit a line, $y = mx + c$, through some noisy data-points:

import numpy as np

x = np.array([0, 1, 2, 3])
y = np.array([-1, 0.2, 0.9, 2.1])

A = np.vstack([x, np.ones(len(x))]).T
A

array([[0., 1.],
[1., 1.],
[2., 1.],
[3., 1.]])

By examining the coefficients, we see that the line should have a gradient of roughly 1 and cut the $y$-axis at, more or less, -1.

We can rewrite the line equation as $y = Ap$, where $A = [[x 1]]$ and $p = [[m], [c]]$. Now use lstsq to solve for $p$:

m, c = np.linalg.lstsq(A, y, rcond=None)[0]
print(m, c)

0.9999999999999999 -0.9499999999999997


Plot the data along with the fitted line:

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 12, 8
plt.rcParams["font.size"] = "17"
plt.plot(x, y, 'o', label='Original data', markersize=10)
plt.plot(x, m*x + c, 'r', label='Fitted line')
plt.legend()
plt.show()

スポンサーリンク

## rcond=Noneの意味¶

np.linalg.lstsq(A, y, rcond=None)のrcond=Noneとは何なのか？こいつが何をしているのかを調べるために、先ずこいつを省いてみた。

m, c = np.linalg.lstsq(A, y)[0]
print(m, c)

0.9999999999999999 -0.9499999999999997

/root/.pyenv/versions/3.7.0/envs/py37/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: rcond parameter will change to the default of machine precision times max(M, N) where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass rcond=None, to keep using the old, explicitly pass rcond=-1.
"""Entry point for launching an IPython kernel.


rcond=Noneは警告をサプレッションするために必要らしい。よく見ると関数の説明書きにも書いてあった。

m, c = np.linalg.lstsq(A, y, rcond=0)[0]
print(m, c)

0.9999999999999999 -0.9499999999999997

スポンサーリンク

## scipy.linalg.lstsq¶

scipy.linalg.lstsqとnumpy.linalg.lstsqを比較する。scipy版についてはこのサイトに以下のように書いてある。

Compute a vector $x$ such that the 2-norm $|b – A x|$ is minimized.

2-ノルム$|b – A x|$が最少になるようなベクトル$x$を算出する。

from scipy.linalg import lstsq

n, d = lstsq(A, y)[0]
print(n, d)

0.9999999999999999 -0.9499999999999997

import matplotlib.pyplot as plt

plt.plot(x, y, 'o', label='Original data', markersize=10)
plt.plot(x, n*x + d, 'r', label='Fitted line')
plt.legend()
plt.show()


スポンサーリンク
スポンサーリンク

フォローする