3.scipy Linalg
title: "3.scipy-linalg" date: 2026-05-24T14:11:09Z
1import pylab as pl
2import numpy as np
3from scipy import linalg
1import matplotlib as mpl
2mpl.rcParams['font.sans-serif'] = ['SimHei']
线性代数-linalg
解线性方程组
1import numpy as np
2from scipy import linalg
3
4m, n = 500, 50
5A = np.random.rand(m, m)
6B = np.random.rand(m, n)
7X1 = linalg.solve(A, B)
8X2 = np.dot(linalg.inv(A), B)
9print (np.allclose(X1, X2))
10%timeit linalg.solve(A, B)
11%timeit np.dot(linalg.inv(A), B)
True
5.38 ms ± 120 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
8.14 ms ± 994 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
1luf = linalg.lu_factor(A)
2X3 = linalg.lu_solve(luf, B)
3np.allclose(X1, X3)
True
1M, N = 1000, 100
2np.random.seed(0)
3A = np.random.rand(M, M)
4B = np.random.rand(M, N)
5Ai = linalg.inv(A)
6luf = linalg.lu_factor(A)
7%timeit linalg.inv(A)
8%timeit np.dot(Ai, B)
9%timeit linalg.lu_factor(A)
10%timeit linalg.lu_solve(luf, B)
50.6 ms ± 1.94 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
3.49 ms ± 306 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
20.1 ms ± 1.42 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
4.49 ms ± 65 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
最小二乘解
1from numpy.lib.stride_tricks import as_strided
2
3
4def make_data(m, n, noise_scale): #❶
5 np.random.seed(42)
6 x = np.random.standard_normal(m)
7 h = np.random.standard_normal(n)
8 y = np.convolve(x, h)
9 yn = y + np.random.standard_normal(len(y)) * noise_scale * np.max(y)
10 return x, yn, h
11
12
13def solve_h(x, y, n): #❷
14 X = as_strided(
15 x, shape=(len(x) - n + 1, n), strides=(x.itemsize, x.itemsize)) #❸
16 Y = y[n - 1:len(x)] #❹
17 h = linalg.lstsq(X, Y) #❺
18 return h[0][::-1] #❻
1x, yn, h = make_data(1000, 100, 0.4)
2H1 = solve_h(x, yn, 120)
3H2 = solve_h(x, yn, 80)
4
5print("Average error of H1:", np.mean(np.abs(h[:100] - h)))
6print("Average error of H2:", np.mean(np.abs(h[:80] - H2)))
Average error of H1: 0.0
Average error of H2: 0.2958422158342371
1#%figonly=实际的系统参数与最小二乘解的比较
2fig, (ax1, ax2) = pl.subplots(2, 1, figsize=(6, 4))
3ax1.plot(h, linewidth=2, label=u"实际的系统参数")
4ax1.plot(H1, linewidth=2, label=u"最小二乘解H1", alpha=0.7)
5ax1.legend(loc="best", ncol=2)
6ax1.set_xlim(0, len(H1))
7
8ax2.plot(h, linewidth=2, label=u"实际的系统参数")
9ax2.plot(H2, linewidth=2, label=u"最小二乘解H2", alpha=0.7)
10ax2.legend(loc="best", ncol=2)
11ax2.set_xlim(0, len(H1));
特征值和特征向量
1A = np.array([[1, -0.3], [-0.1, 0.9]])
2evalues, evectors = linalg.eig(A)
3
4print(evalues)
5print(evectors)
[1.13027756+0.j 0.76972244+0.j]
[[ 0.91724574 0.79325185]
[-0.3983218 0.60889368]]
1#%figonly=线性变换将蓝色箭头变换为红色箭头
2points = np.array([[0, 1.0], [1.0, 0], [1, 1]])
3
4def draw_arrows(points, **kw):
5 props = dict(color="blue", arrowstyle="->")
6 props.update(kw)
7 for x, y in points:
8 pl.annotate("",
9 xy=(x, y), xycoords='data',
10 xytext=(0, 0), textcoords='data',
11 arrowprops=props)
12
13draw_arrows(points)
14draw_arrows(np.dot(A, points.T).T, color="red")
15draw_arrows(evectors.T, alpha=0.7, linewidth=2)
16draw_arrows(np.dot(A, evectors).T, color="red", alpha=0.7, linewidth=2)
17
18ax = pl.gca()
19ax.set_aspect("equal")
20ax.set_xlim(-0.5, 1.1)
21ax.set_ylim(-0.5, 1.1);
1np.random.seed(42)
2t = np.random.uniform(0, 2*np.pi, 60)
3
4alpha = 0.4
5a = 0.5
6b = 1.0
7x = 1.0 + a*np.cos(t)*np.cos(alpha) - b*np.sin(t)*np.sin(alpha)
8y = 1.0 + a*np.cos(t)*np.sin(alpha) - b*np.sin(t)*np.cos(alpha)
9x += np.random.normal(0, 0.05, size=len(x))
10y += np.random.normal(0, 0.05, size=len(y))
1D = np.c_[x**2, x*y, y**2, x, y, np.ones_like(x)]
2A = np.dot(D.T, D)
3C = np.zeros((6, 6))
4C[[0, 1, 2], [2, 1, 0]] = 2, -1, 2
5evalues, evectors = linalg.eig(A, C) #❶
6evectors = np.real(evectors)
7err = np.mean(np.dot(D, evectors)**2, 0) #❷
8p = evectors[:, np.argmin(err) ] #❸
9print (p)
[-0.55214278 0.5580915 -0.23809922 0.54584559 -0.08350449 -0.14852803]
1#%figonly=用广义特征向量计算的拟合椭圆
2def ellipse(p, x, y):
3 a, b, c, d, e, f = p
4 return a*x**2 + b*x*y + c*y**2 + d*x + e*y + f
5
6X, Y = np.mgrid[0:2:100j, 0:2:100j]
7Z = ellipse(p, X, Y)
8pl.plot(x, y, "ro", alpha=0.5)
9pl.contour(X, Y, Z, levels=[0]);
奇异值分解-SVD
1r, g, b = np.rollaxis(pl.imread("vinci_target.png"), 2).astype(float)
2img = 0.2989 * r + 0.5870 * g + 0.1140 * b
3img.shape
(505, 375)
1U, s, Vh = linalg.svd(img)
2print(U.shape)
3print(s.shape)
4print(Vh.shape)
(505, 505)
(375,)
(375, 375)
1#%fig=按从大到小排列的奇异值
2pl.semilogy(s, lw=3);
C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\mathtext.py:854: MathTextWarning: Font 'default' does not have a glyph for '-' [U+2212]
MathTextWarning)
C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\mathtext.py:855: MathTextWarning: Substituting with a dummy symbol.
warn("Substituting with a dummy symbol.", MathTextWarning)
1def composite(U, s, Vh, n):
2 return np.dot(U[:, :n], s[:n, np.newaxis] * Vh[:n, :])
3
4print (np.allclose(img, composite(U, s, Vh, len(s))))
True
1#%fig=原始图像、使用10、20、50个向量合成的图像（从左到右）
2img10 = composite(U, s, Vh, 10)
3img20 = composite(U, s, Vh, 20)
4img50 = composite(U, s, Vh, 50)
1%array_image img; img10; img20; img50
UsageError: Line magic function `%array_image` not found.
1pl.imshow(img)
<matplotlib.image.AxesImage at 0x16bf1557b00>
1pl.imshow(img10)
<matplotlib.image.AxesImage at 0x16bf186fc88>
1pl.imshow(img20)
<matplotlib.image.AxesImage at 0x16bf1d99cc0>
1pl.imshow(img50)
<matplotlib.image.AxesImage at 0x16bf1de9da0>
