理解特异值分解(SVD)

svd_intro.py

本习题的程序文件

在大学的线性代数课程中通常不会涉及特异值分解(SVD)相关的内容,然而SVD却是线性代数中一个十分重要的内容。在本习题中我们通过一系列的练习帮助读者理解SVD。

有关SVD的详细介绍,请参考下面的链接:

http://www.uwlax.edu/faculty/will/svd/

Introduction to the Singular Value Decomposition

坐标变换

我们可以用一个二维矢量表示二维平面上的某个点,而一个2x2的矩阵则可以用来表示二维平面上的某种线性坐标变换。当用变换矩阵和矢量相乘,我们就可以得到变换之后的矢量。为了演示坐标变换的效果,我们首先编写一个用matplotlib绘制变换前后坐标点的函数。

def plot_arrows(p0, p1):
    axe = pl.gca()
    axe.set_aspect("equal")

    for i in xrange(p0.shape[1]):
        from_ = p0[:, i]
        delta = p1[:, i] - from_
        if delta[0]!=0 and delta[1]!=0:
            axe.arrow(from_[0], from_[1], delta[0], delta[1], head_width=0.1, length_includes_head=True)
    axe.scatter(p0[0,:], p0[1,:], s=50, c="blue")
    axe.scatter(p1[0,:], p1[1,:], s=50, c="red")

plot_arrows(p0, p1)绘制p0和p1所表示的点,以及从p0到p1的箭头。其中p0和p1是形状为(2, N)的数组,N为点的个数。请编写程序用plot_arrows()绘制矩阵m对单位圆上的点的变换结果。其效果如【图:矩阵m对单位圆的变换】所示。

m = np.array([
    [ 1.0, 3.0],
    [-3.0, 2.0]])
【你的程序】计算size列表,它的每个小标都保存增长到此下标时size列表的大小
/tech/static/books/scipyex/_images//svd_intro_01.png

矩阵m对单位圆的变换

标准正交基(二维)

在n维空间中,n个互相垂直且长度为1的向量构成此空间中的标准正交基。我们可以把这n个向量想象成n维空间的一组垂直坐标系。对于二维空间来说,将X轴和Y轴的单位向量绕原点旋转\theta即可得到一组标准正交基。即:

v_1 = \begin{bmatrix}
\cos \theta\\
\sin \theta\\
\end{bmatrix},  v_2 = \begin{bmatrix}
-\sin \theta \\
\cos \theta \\
\end{bmatrix}

当使用以标准正交基为行向量的矩阵进行坐标变换时,它能将标准正交基变换成单位向量。因为此矩阵将正交基向量变换成单位向量,因此它也被称为对准矩阵。在【图:对准矩阵的变换演示】中,蓝色箭头构成单位向量,而红色箭头为一组标准正交基。通过标准正交基构成的对准矩阵,可以把红色点所在的坐标变换为蓝色点所在的坐标。也就是当红色坐标系旋转到蓝色坐标系时,红色点也正好和蓝色点重合。也可以将这种变换理解为将X-Y坐标系中的坐标变换为v_1- v_2坐标系中的坐标。图中,红色点在蓝色坐标系中的坐标在经过变换之后,将变为蓝色点在蓝色坐标系中的坐标,它等于红色点在红色坐标系中的坐标。

/tech/static/books/scipyex/_images//svd_intro_02.png

对准矩阵的变换演示

下面是绘制【图:对准矩阵的变换演示】的程序。

def make_perpframes2d(theta):
    v1 = np.array([ np.cos(theta), np.sin(theta)])
    v2 = np.array([-np.sin(theta), np.cos(theta)])
    return v1, v2

def plot_axis(axis, name, color):
    axe = pl.gca()
    axe.arrow(0, 0, axis[0], axis[1], head_width=0.05,
              length_includes_head=True, color=color)
    axe.text(axis[0], axis[1], name, fontsize=16)

v1, v2 = make_perpframes2d(np.pi/6)
m = np.vstack((v1, v2))
p0 = np.array([[0.3,0.4]])
p1 = np.dot(m, p0.T).T

pl.figure()
pl.subplot(111)
axe = pl.gca()
axe.set_aspect("equal")
plot_axis((1,0), "$x$", "blue")
plot_axis((0,1), "$y$", "blue")
plot_axis(v1, "$v_1$", "red")
plot_axis(v2, "$v_2$", "red")
pl.scatter(p0[:,0], p0[:, 1], c="red")
pl.scatter(p1[:,0], p1[:, 1], c="blue")
axe.set_xlim(-0.8,1.2)
axe.set_ylim(-0.2,1.2)

当使用以标准正交基为列变量的矩阵进行变换时,它的变换效果正好和对准矩阵相反,通常被称为悬挂矩阵。因为它将单位向量“悬挂”到标准正交基上。

请仿照对准矩阵的变换示例,编写悬挂矩阵变换的演示程序。

标准正交基(三维)

二维平面上的标准正交基只有一个自由度,只要给定一个旋转角度参数,即可唯一决定一组标准正交基。这是因为两个向量v_1 = \begin{bmatrix}a \\ b \end{bmatrix}v_2 = \begin{bmatrix}c \\ d \end{bmatrix}虽然有四个未知数,但是根据正交基的定义,它们必须满足以下三个条件:

a^2 + b^2 = 1 \\
c^2 + d^2 = 1 \\
a \cdot c + b \cdot d = 0

前两个方程决定了两个向量的长度为1,而最后一个方程决定了两个向量是正交的。

对于三维空间中的标准正交基,有9个未知数,它们必须满足三个长度为1的方程,和三个正交方程,因此三维空间中的标准正交基有三个自由度。

请编写用VPython绘制演示三维对准矩阵变换的程序。程序的输出应与【图:三维空间中的对准矩阵变换演示】类似。

/tech/static/books/scipyex/_images//svd_intro_03.png

三维空间中的对准矩阵变换演示

svd_intro_3d_transform.py

用VPython演示三维空间的标准正交基和对准矩阵变换

import pylab as pl
import numpy as np
from scipy.optimize import fsolve

def make_perpframes3d(a,b,d):
    """通过a,b,d三个值,计算一组三维空间的标准正交基:
    v1 = (a,b,c)
    v2 = (d,e,f)
    v3 = (k,i,j)
    其中a,b,d为已知数, c,e,f,k,i,j为未知数。
    """
    【你的程序】

def show_preframes(m, color, axis_label):
    for l, v in zip(axis_label, m):
        arrow(pos=(0,0,0), axis=v, color=color, shaftwidth=0.02)
        label(text=l, pos=v*1.1, opacity=0, box=False, color=(0,0,0), height=20)

if __name__ == "__main__":
    from visual import *
    scene.background = 1,1,1
    m = make_perpframes3d(0.2, 0.5, -0.3)
    p0 = np.array([-0.5, 0.6, 0.7])
    p1 = np.dot(m, p0.T).T
    show_preframes(np.identity(3), (0, 1, 0), ("x","y","z"))
    show_preframes(m, (1, 0, 0), ("v1","v2","v3"))
    sphere(pos=p0,radius=0.03, color=(1,0,0))
    sphere(pos=p1,radius=0.03, color=(0,1,0))
可以使用scipy.optimize.fsolve()对标准正交基的方程组求解。这样可以通过指定的三个自由度的值计算一组标准正交基。

內容目录

上一个主题

PID控制

本页

loading...