numpy积累

这篇文章主要解决的问题是，Numpy 提供了许多线性代数的工具，但是我因为线性代数学的不好，因此总是不会用。

基本公式 of linalg

numpy.linalg.svd

函数定义：

linalg.svd(a, full_matrices=True, compute_uv=True, hermitian=False)

Singular Value Decomposition.

When _a_ is a 2D array, and full_matrices=False, then it is factorized as u @ np.diag(s) @ vh = (u * s) @ vh, where _u_ and the Hermitian transpose of _vh_ are 2D arrays with orthonormal columns and _s_ is a 1D array of _a_’s singular values. When _a_ is higher-dimensional, SVD is applied in stacked mode as explained below.

numpy.linalg.svd — NumPy v1.24 Manual

SVD 的作用是将 2D 矩阵进行分解。对于高维的矩阵，之后再研究。矩阵 $A$，可以写成 SVD 分解形式：$A=USV^H$，其中 $A=a,U=u,S=np.diag(s), V^H=vh$，其中 $s$ 是矩阵 $A$ 的所有奇异值 $\sigma_i$，有 $\sigma_i=\sqrt{\lambda_i}$，$vh$ 的行是 $A^HA$ 的特征向量，$u$ 的列是 $AA^H$ 的特征向量，而且 $u,vh$ 都是酉矩阵，转置相乘为单位矩阵。这里的上标的意思是转置。

特征值分解，奇异值分解（SVD） - 知乎

注意得到的特征向量都不一定是 $A$ 的。因此 DLT 算法中：

我们想求的是 $A^\top A$ 的特征向量，然后就有 $h^\top A^\top Ah=\sigma_i h^\top h=\sigma_i |h|^2 = \sigma_i$，也就是说要使得特征值最小，那不就取最小特征值的特征向量，也就是 $vh[-1]$。

如果 full_matrices=True，则三个返回值的维度为 mxm, min(m,n), nxn，这时候想要重建原本的矩阵，需要将中间的矩阵还原：

u, s, vh = np.linalg.svd(a, full_matrices=True)
u.shape, s.shape, vh.shape # ((9, 9), (6,), (6, 6))
np.allclose(a, np.dot(u[:, :6] * s, vh)) # True
smat = np.zeros((9, 6), dtype=complex)
smat[:6, :6] = np.diag(s)
np.allclose(a, np.dot(u, np.dot(smat, vh)))# True

如果使用 full_matrices=False，则可以直接相乘。

u, s, vh = np.linalg.svd(a, full_matrices=False)
u.shape, s.shape, vh.shape # ((9, 6), (6,), (6, 6))
np.allclose(a, np.dot(u * s, vh)) # True
smat = np.diag(s)
np.allclose(a, np.dot(u, np.dot(smat, vh))) # True

numpy.linalg.pinv

numpy.linalg.pinv — NumPy v1.24 Manual

linalg.pinv(a, rcond=1e-15, hermitian=False)

计算矩阵的伪逆。伪逆是一个没有学过的概念，参考 Moore–Penrose inverse - Wikipedia。只需要记住最简单的定义：

对于一个有解的方程 $Ax=b$，伪逆可以给出一个矩阵 $A^+$，使得一个解 $\bar{x}=A^+b$ 成立。
事实上，如果对 $A$ 进行分解，$A=Q_1\Sigma Q_2^\top$，可以认为 $A^+=Q_2\Sigma^+Q_1^\top$，其中 $\Sigma^+$ 是 $\Sigma$ 的各特征值的倒数。

因此，可以用这个函数进行问题的求解，比如 Photometric Stereo 问题：

求解出来 $\tilde{n}$ 即可。

numpy.linalg.norm

numpy.linalg.norm — NumPy v1.24 Manual
作用是对选中的轴进行归一化，得到范数，返回值是结果。不是原地变化。

linalg.norm(x, ord=None, axis=None, keepdims=False)

如果不选择 axis，则对所有的元素求范数。
ord 的作用如下：
For values of ord < 1, the result is, strictly speaking, not a mathematical ‘norm’, but it may still be useful for various numerical purposes.
The following norms can be calculated:

ord	norm for matrices	norm for vectors
None	Frobenius norm	2-norm
fro	Frobenius norm	–
nuc	nuclear norm	–
inf	max(sum(abs(x), axis=1))	max(abs(x))
-inf	min(sum(abs(x), axis=1))	min(abs(x))
0	–	sum(x != 0)
1	max(sum(abs(x), axis=0))	as below
-1	min(sum(abs(x), axis=0))	as below
2	2-norm (largest sing. value)	as below
-2	smallest singular value	as below
other	–	sum(abs(x)ord) (1./ord)

The Frobenius norm is given by:

$$\newcommand{\d}{\mathrm{d}} \newcommand{\bs}[1]{\boldsymbol{ #1}} ||A||_F=[\sum_{i,j}\text{abs}(a_{i,j})^2]^{1/2}$$

The nuclear norm is the sum of the singular values.
Both the Frobenius and nuclear norm orders are only defined for matrices and raise a ValueError when x.ndim != 2.