偏题: 矩阵的逆
注: 这部分并不是课上介绍的, 实际上也不应该自己来写, 对于 Python 来说, 使用 numpy
调包, 用 np.inv()
的方法是最快的.
因为我 吃饱了撑的, 所以有了下面的矩阵的逆的部分.
不过因为与课程没啥关系, 所以用的是个人比较喜欢的 Common Lisp.
Gauss 消元法
最简单的版本就是使用高斯消元法, 实现的方法看线代即可, 不过这里并不详细介绍该方法, 因为实现起来有点麻烦.
LU 分解求逆矩阵
思路:
- 将矩阵分解为 $A = L U$, 其中 $L$ 为下三角矩阵, 而 $U$ 为上三角矩阵
- 对于三角矩阵的逆的求解和 Gauss 消元法的最后的解方程的操作一致
LU 分解
算法如下图示:
u u u u u u u u u u u u u u u
l u u u u l u u u u l u u u u
l l # r r <- k => l l * * * => l l * * *
l l r r r l l r r r l l % r r
l l r r r (2) l l r r r (3) l l % r r (4)
以下假设矩阵的脚标从零计数, 假设矩阵的大小为 rows * cols
, 并且是在矩阵原位上进行计算, 即最终的 LU 矩阵是用一个矩阵来进行记录的.
k = 0
- 计算第
k
层的 u
和 l
, 首先得到 $a_{k,k}$
- 计算上三角矩阵 U 的
k
层: $a{i,k} = a{i,k} / a_{k,k}, i \in (k, rows)$
- 计算下三角矩阵 L 的
k
层和剩余的 R 的部分: $a{i,j} = a{i,j} - a{i,k} \times a{k,j}$
这里使用 Common Lisp 代码为例:
(defun matrix-lu-decomp (matrix &optional (compress? t))
"LU decomposition of list `matrix'.
Return values are L-matrix and U-matrix.
If `compress?' is non-nil, return LU in a single matrix,
otherwise, return values are L, U matrix seperately.
Example:
(matrix-lu-decomp '((4 3)
(6 3)))
;; => ((4 3)
(1.5 -1.5))
Algorithm:
u u u u u u u u u u u u u u u
l u u u u l u u u u l u u u u
l l # r r <- k => l l * * * => l l * * *
l l r r r l l r r r l l % r r
l l r r r (a) l l r r r (b) l l % r r (c)
"
(multiple-value-bind (rows cols)
(matrix-size matrix)
(let ((s (min rows cols))
(matrix (copy-tree matrix)))
(macrolet ((a (i j)
`(nth ,j (nth ,i matrix))))
(loop for k below s
for x = (/ 1.0 (a k k)) ;; 1/# in (a)
do (loop for i from (1+ k) below rows ;; * in (b)
do (setf (a i k) (* x (a i k)))) ;; a_ik <- a_ik / a_kk
do (loop for i from (1+ k) below rows ;; % in (c)
do (loop for j from (1+ k) below cols
do (setf (a i j) ;; a_ij <- a_ij - a_ik * a_kj
(- (a i j) (* (a i k) (a k j))))))
finally
(return ;; if `compress?' return the compressed matrix
(if compress? matrix (decompress-lu-matrix matrix))))))))
三角矩阵的逆的求解
以常见的上三角为例, 假设大小为 n
:
| * * * * | x1 | | 0 |
| 0 * * * | x2 | = | . | => X_i = ( ..... 1/a_{ii} 0 ....... 0)
| 0 0 * * | .. | | 1 | <- i \ i-1 / i \ n-i-1 /
| 0 0 0 * | xn | | 0 |
只需要求解方程 $U X_i = (0, 0, \dots, 1_i, \dots, 0)$, 然后把 $X_i$ 拼合成一个矩阵. 拼合的过程如下:
i = 0
vec = [1/a_{i,i}, 0, ..., 0]
, 0
有 n - i - 1
个
j = i - 1
, weight = [a_{k,j}, k = (j+1)...n]
vec = [- dot(weight, vec) / a_{j,j}] + vec
这里使用 Common Lisp 的代码为例:
(defun u-inverse (lu-matrix)
"Calculate the inverse matrix of upper triangle list matrix on `lu-matrix'.
Return value are the inverse matrix of `matrix'.
Algorithm:
| * * * * | x1 | | 0 |
| 0 * * * | x2 | = | . | => X_i = ( ..... 1/a_{ii} 0 ....... 0)
| 0 0 * * | .. | | 1 | <- i \ i-1 / i \ n-i-1 /
| 0 0 0 * | xn | | 0 |
"
(let ((size (length lu-matrix)))
(loop for i below size
;; calculate the x_i for U X = I
for vec = (cons (/ 1 (u-matrix-at lu-matrix i i))
(make-list (- size i 1) :initial-element 0))
;; calculate the x_k, k < i
do (loop for j from (1- i) downto 0
for weight = (loop for k from (1+ j) below size
collect (u-matrix-at lu-matrix k j))
for dot = (loop with dot = 0
for k from j below i
do (incf dot (u-matrix-at lu-matrix k j))
finally (return dot))
do (setf vec (cons (/ (- dot) (u-matrix-at lu-matrix j j)) vec)))
collect vec into mat
finally (return (matrix-transpose mat)))))
(defun l-inverse (lu-matrix)
"Calculate the inverse matrix of lower triangle list matrix `matrix'.
Return value are the inverse matrix of `matrix'.
Algorithm:
| * 0 0 0 | x1 | | 0 |
| * * 0 0 | x2 | = | . | => X_i = ( 0 ..... 0 1/a_{ii} ... )
| * * * 0 | xj | | 1 | <- i \ i-1 /
| * * * * | xn | | . |
^k
"
(let ((size (length lu-matrix)))
(loop
;; calculate the X_i as `reverse-vec'
for i below size
;; first i-1 element is 0
;; the `reverse-vec' is for X_i, but reversed
for reverse-vec = (make-list i :initial-element 0)
;; i element is 1/a_ii
do (push (/ 1.0 (l-matrix-at lu-matrix i i)) reverse-vec)
;; calculate the x_j, j > i in X_i
;; x_j = - (a_{k,j} * x_k) / a_{j,j}, k < j
do (loop for j from (1+ i) below size
for dot = (loop with dot = 0
for k from (1- j) downto i
for x_k in reverse-vec
do (incf dot (* (l-matrix-at lu-matrix k j) x_k))
finally (return dot))
do (push (- (/ dot (l-matrix-at lu-matrix j j)))
reverse-vec))
collect (nreverse reverse-vec) into mat
finally (return (matrix-transpose mat)))))
于是可以计算逆矩阵:
(defun matrix-inverse (matrix)
"Return the inverse matrix for list matrix `matrix'.
Note that there maybe precision loss while in computation.
Algorithm:
M = L U => M^{-1} = U^{-1} L^{-1}
"
(let ((lu-matrix (matrix-lu-decomp matrix)))
(mat-mat (u-inverse lu-matrix)
(l-inverse lu-matrix))))
BLAS/LAPACK
虽然但是, 直接这么做还是有点性能低下呢. 虽然也不是不能用多线程之类的方式进行加速计算, 但是如果能够调库的话.
见 Fortran PDGETRF, DGETRF 代码, 或是直接用其他语言的 binding 吧.
顺带一提, 如果是使用 CFFI 进行调用的话, 只需要用 dgetrf_
这个函数即可.
以 Common Lisp 中的 CFFI 为例:
(CFFI:DEFCFUN ("dgetrf_" %DGETRF) :VOID (M FORTRAN-INT) (N FORTRAN-INT)
(A CFFI-FNV-DOUBLE) (LDA FORTRAN-INT) (IPIV CFFI-FNV-INT32)
(INFO FORTRAN-INT))
(EXPORT '%DGETRF :LAPACK-CFFI)
(代码来源于 cl-blapack)