线性代数学习笔记二

1 线性方程组

矩阵记号是为解方程组带来方便。解方程组，消元法。

三种基本变换对应于增广矩阵的下列变换:

行初等变换

(倍加变换 replacement) 把某一行换成它本身与另一行的倍数的和
(对换变换 interchange) 把两行对换
(倍乘变换 scaling) 把某一行的所有元素乘以同一个非零数

行变换可应用于任何矩阵.如果一个矩阵可以经过一系列行初等变换变成另一个矩阵，则称这两个矩阵是行等价的。

行变换是可逆的。

线性方程组的两个基本问题：

方程组是否相容，即它是否至少有一个解？
若它有解，是否只有一个解，即解是否唯一？

1.1 行化简与阶梯型矩阵

非零行或列指矩阵中至少包含一个非零元素的行或列，非零行的先导元素是指该行中最左边的非零元素。

一个矩阵称为阶梯型，则有以下三个性值: 1 每一非零行在每一零行之上 2 某一行的先导元素所在的列位于前一行先导元素的右面 3 某一先导元素所在列下方元素都是零。若一个阶梯型矩阵还满足以下性质，则称它为简化阶梯形： 4 每一非零行的先导元素是1 5 每一先导元素1是该元素所在列的唯一非零元素

一个矩阵可以行化简变为阶梯形矩阵，但不同的方法可以化为不同的阶梯形矩阵。一个矩阵只能化为唯一的简化阶梯形矩阵。

每个矩阵行等价于唯一的简化阶梯形矩阵。大部分矩阵程序用RREF作为简化阶梯形的缩写，有些用REF作为阶梯形的缩写。

1.1.1 主元位置

矩阵经过行变换化为阶梯形后，经进一步的行变换将矩阵化为简化阶梯形时，先导元素的位置并不改变。

矩阵中的主元位置是A中对应于它的阶梯形中先导元素的位置。主元列是A的含有主元位置的列。

用不同的行变换产生不同的主元。

行简化算法，向前步骤，向后步骤。

1.1.2 线性方程组的解

行简化算法应用于方程组的增广矩阵时，可以得出线性方程组解集的一种显式示法。

对应主元列的变量称为基本变量，其它变量称为自由变量。自由变量是指它可取任意值。

解集的参数表示

有自由变量的方程有无穷多个解

线性方程组相容的充要条件是增广矩阵的最右列不是主元列。就是说，增广矩阵的阶梯形没有 [0 … 0 b] b ≠ 0的行，如果线性方程组相容，它的解集可能有两种情形:

当没有自由变量时，有唯一解。
若至少有一个自由变量，有无穷多解。

#!/usr/bin/env plk

;;如果使用lein exec 需要添加alembic依赖,jvm clojure启动太慢了,
;;改用planck, clojurescript REPL,速度快

;(require 'alembic.still)
;(alembic.still/distill '[org.clojure/math.numeric-tower "0.0.4"])
;(require '[clojure.math.numeric-tower :as math])

;;; matrix
(defn make-matrix
  "构造一个矩阵,
  array为元素内容
  len为每行有几个元素"
  [array len]
  (partition len array))

(defn cols-of-matrix
  "矩阵的列数"
  [m]
  (count (first m)))

(defn rows-of-matrix
  "矩阵的行数"
  [m]
  (count m))

(defn drop-ncol-matrix
  "丢弃矩阵m的前n列"
  [m n]
  (map #(drop n %) m))

(defn print-matrix
  "打印一个矩阵"
  [m]
  (println "--rows:" (rows-of-matrix m) ",cols:" (cols-of-matrix m) "--" )
  (doseq [row m]
    (print "|")
    (doseq [v row]
      (print v "	" ))
    (println "|")))

(defn transpose-matrix
  "矩阵转置，就是行变为列，列变为行"
  [m]
  (apply mapv vector m))

(defn add-matrix
  "加两个矩阵"
  [m1 m2]
  (assert (and
           (= (cols-of-matrix m1) (cols-of-matrix m2))
           (= (rows-of-matrix m1) (rows-of-matrix m2)))
          "matrix length not match")
  (map #(map + %1 %2) m1 m2))

(defn scal-matrix
  "放大矩阵m"
  [scal-n matrix]
  (map (fn [row]
         (map #(* scal-n %1) row)) matrix))

(defn dot-product
  "两个向量的点积"
  [v1 v2]
  (reduce + (map * v1 v2)))

(defn mul-matrix
  "乘2个矩阵"
  [m1 m2]
  (assert (= (cols-of-matrix m1) (rows-of-matrix m2))
          "matrix not match, can't mul")
  (let [tm2 (transpose-matrix m2)]
    (map (fn [row-m1]
           (map (fn [col-m2]
                  (dot-product row-m1 col-m2))
                tm2))
         m1)))

(defn col-of-matrix
  "返回一个矩阵m的第c列"
  [m c]
  (map #(nth % c) m))

(defn row-of-matrix
  "返回一个矩阵m的第r行"
  [m r]
  (nth m r))

(defn val-of-matrix
  "返回矩阵m位置(row,col)的值"
  [m row col]
  (nth (nth m row) col))

(defn pivot-of-vec
  "返回一个向量中主元的索引及主元值"
  [v]
  (let [[idx _](apply max-key second (map-indexed #(vector %1 (Math/abs %2)) v))]
    [idx (nth v idx)]))

(defn replacement-row
  "倍加变换"
  [m r1 r2 n]
  (println "add " n " times row " r2 " to row " r1 )
  (let [row2 (row-of-matrix m r2)
        scal-row2 (map #(* n %1) row2)]
    (map-indexed (fn [idx row]
                   (if (== idx r1)
                     (map + row scal-row2)
                     row))
                 m)))

(defn interchange-row
  "交换矩阵m的行r1和r2"
  [m r1 r2]
  (let [row1 (row-of-matrix m r1)
        row2 (row-of-matrix m r2)]
    (println "interchange row " r1 " and " r2)
    (map-indexed (fn [idx r]
                   (cond
                     (= idx r1) row2
                     (= idx r2) row1
                     :else r)) m)))

(defn scaling-row
  "放大矩阵m第r行n倍,倍乘r行"
  [m r n]
  (println "scaling row " r " with " n)
  (map-indexed (fn [idx row]
                 (if (== idx r)
                   (map #(* n %1) row)
                   row)) m))

(defn echelon-form
  "对矩阵m行化简为阶梯形矩阵"
  [m]
  (let [cols (cols-of-matrix m)
        rows (rows-of-matrix m)]
    (loop [row-idx 0
           col-idx 0
           result m]
      (println "loop forward phase row" row-idx " col-idx" col-idx "matrix:")
      (print-matrix result)
      (if (= row-idx (- rows 1))
        result
        (let [rest (take-last (- rows row-idx) result)
              [pivot-idx pivot col-idx]
                                        ;找到主元列
              (loop [idx col-idx]
                (let [[pivot-idx pivot] (pivot-of-vec (col-of-matrix rest idx))]
                  (if (== 0 pivot)
                    (recur (inc idx))
                                        ;这里的pivot-idx加上row-idx是因为此矩阵是截取过的
                    [(+ row-idx pivot-idx) pivot idx])))
;              a (println "pivot of col" col-idx " pivot-row" pivot-idx " pivot " pivot )
                                        ;是否交换行

              result (if-not (= row-idx pivot-idx)
                       (interchange-row result row-idx pivot-idx)
                       result)
                                        ;倍加消元
              result (reduce (fn [r idx]
                               (let [rv (val-of-matrix r idx col-idx)]
                                 (if (== 0 rv)
                                   r
                                   (replacement-row r idx row-idx (/ (- rv) pivot))))
                               ) result (range (inc row-idx) rows))]
          (recur (inc row-idx)
                 (inc col-idx)
                 result))))))

(defn pivot-of-row
  "获得矩阵某行的主元,返回[c n]
  c为主元所在列
  n为主元值
  "
  [m r]
  (->> (row-of-matrix m r)
       (map-indexed vector)
       (filter (fn [[idx v]] (not (== v 0))))
       first))

(defn reduced-echelon-form
  "简化阶梯形矩阵"
  [m]
  (loop [row-idx (- (rows-of-matrix m) 1)
         result (echelon-form m)]
    (println "loop backward phase,row " row-idx )
    (print-matrix result)
    (if (< row-idx 0)
      result
      (let [[pivot-col pivot] (pivot-of-row result row-idx)
            _ (println "pivot at (" row-idx "," pivot-col ") = " pivot)
            result (if (nil? pivot)
                     result
                     (loop [result (if-not (== 1 pivot)
                                     (scaling-row result row-idx (/ 1 pivot))
                                     result)
                            r2-idx (dec row-idx)]
                                        ;倍加消元
                       (if (< r2-idx 0)
                         result
                         (let [v (val-of-matrix result r2-idx pivot-col)]
                           (if (== 0 v)
                             (recur result (dec r2-idx))
                             (recur (replacement-row result r2-idx row-idx (- v))
                                    (dec r2-idx)))))))]
        (recur (dec row-idx) result)))))

(def a (make-matrix [3 -1 2 0] 2))
(def b (make-matrix [-7 2 3 5] 2))
(def c (make-matrix [-7 2 3 5 6 8 9 2 3 4] 5))
(print-matrix a)
(print-matrix b)
(print-matrix c)
(echelon-form c)
(reduced-echelon-form c)
(print-matrix (add-matrix a b))
(print-matrix (add-matrix a (scal-matrix -1 b)))
(print-matrix (add-matrix
               (make-matrix [1 2 3 4 5 6 7 8 9] 3)
               (make-matrix [-10 -100 -1000 1 0 0 1 0 1] 3)))
(print-matrix (mul-matrix
               (make-matrix [2 -3 7 5] 2)
               (make-matrix [10 -8 12 -2] 2)))
(reduced-echelon-form (make-matrix [1 -2 1 0
                                    0 2 -8 8
                                    5 0 -5 10] 4))
;;; AB <> BA
(print-matrix (mul-matrix
               (make-matrix [1 2 3 4] 2)
               (make-matrix [5 6 7 8] 2)))
(print-matrix (mul-matrix
               (make-matrix [5 6 7 8] 2)
               (make-matrix [1 2 3 4] 2)))

(print-matrix (mul-matrix
               (make-matrix [3 1 2 -2 0 5] 3)
               (make-matrix [-1 3 0 5 1 5] 2)))

(def a (make-matrix [1 -2 1 0 0 2 -8 8 5 0 -5 10] 4))
(print-matrix a)
(print-matrix (replacement-row a 2 0 (/ -5 1)))
(print-matrix (scaling-row a 2 (/ 1 5)))
(print-matrix (interchange-row a 0 2))

(def a (make-matrix [0 3 -6 6 4 -5
                     3 -7 8 -5 8 9
                     3 -9 12 -9 6 15] 6))
(print-matrix a)
(print-matrix (echelon-form a))
(reduced-echelon-form a)
(reduced-echelon-form (make-matrix [1 0 -5 1
                                    0 1 1 4
                                    0 0 0 0 ] 4))
(def b (make-matrix [1 6 2 -5 -2 -4
                     0 0 2 -8 -1 3
                     0 0 0 0 1 7] 6))
(reduced-echelon-form b)

(reduced-echelon-form (make-matrix [3 -9 12 -9 6 15
                                    0 2 -4 4 2 -6
                                    0 0 0 0 1 4] 6))
(reduced-echelon-form (make-matrix [1 -3 -5 0
                                    0 1 1 3] 4))
(defn vec-i
  "返回一个向量,
  长度为n,默认值为0,
  第i个元素值为m,"
  [m i n]
  (-> (repeat n 0)
      (vec)
      (assoc-in [i] m)))

;;单位矩阵I  IA=A AI=A 方阵
(defn identity-matrix
  "返回nxn的单位矩阵"
  [n]
  (vec (map-indexed #(vec-i %2 %1 n) (repeat n 1))))

(defn augment-matrix
  "增广矩阵"
  [m1 m2]
  (assert (= (rows-of-matrix m1)
             (rows-of-matrix m2))
          "augment matrix rows mismatch")
  (map concat m1 m2))

(augment-matrix (make-matrix [1 2 3 4 5 6] 3)
                (make-matrix [7 8 9 10 11 12] 3)
                )

(defn invert-matrix
  "获得矩阵的逆矩阵"
  [m]
  (let [rn (rows-of-matrix m) ]
    (assert (= rn (cols-of-matrix m))
            "invert not a square matrix")
    (-> (reduced-echelon-form (augment-matrix m (identity-matrix rn)))
        (drop-ncol-matrix rn))))

(print-matrix (identity-matrix 4))

;; A'A = I ; A'是A的逆
(print-matrix (invert-matrix (make-matrix [2 5 -3 -7] 2)))
(print-matrix (invert-matrix (make-matrix [0 1 2
                                           1 0 3
                                           4 -3 8] 3)))

(defn add-vec
  "加两个向量"
  [v1 v2]
  (map + v1 v2))

(defn scal-vec
  "标量乘法"
  [n v]
  (map  #(* n  %1) v))

(add-vec [4 -8] [-6 15])
(scal-vec 4 [4 -8])
(add-vec (scal-vec 4 [1 -2]) (scal-vec -3 [2 -5]))
(reduced-echelon-form (make-matrix [1 2 7 -2 5 4 -5 6 -3] 3))
(reduced-echelon-form (make-matrix [-5 6 -3 1 2 7 -2 5 4 ] 3))
(add-vec (scal-vec 100 [0.45 0.25 0.15])
         (scal-vec 20 [0.4 0.3 0.15]))
(def v1 [20 550])
(def v2 [30 500])
(reduced-echelon-form (make-matrix [20 30 150
                                    550 500 2825] 3))
(reduced-echelon-form (make-matrix [27.6 30.2 162
                                    3100 6400 23610
                                    250 360 1623] 3))
(mul-matrix (make-matrix [1 2 -1 0 -5 3] 3)
            (make-matrix [4 3 7] 1))

1.2 向量方程

仅含一列的矩阵称为列向量，或简称向量。

所有两个元素的向量的集记为R², R表示向量中的元素是实数，指数2表示每个向量包含2个元素。

给定Rⁿ 中向量v₁, v₂, …, v_p 和标量 c₁, c₂, …, c_p ,向量 [ y = c_1v_1+...+c_pv_p ] 称为向量v₁, v₂, …, v_p 以c₁, c₂, …, c_p 为权的线性组合。

向量方程 [ x_1a_1+x_2a_2+...+x_na_n=b] 和增广矩阵为

egin{equation} label{eq:m1} [a_1 quad a_2 quad ... quad a_n quad b] end{equation}

的线性方程组有相同的解集。特别地， b 可以表示为 a₁, a₂, …, a_n 的线性组合，当且仅当对应于 eqref{eq:m1} 式的方程组有解.

1.3 矩阵方程

若A是mxn矩阵，它的各列为a₁, …, a_n. 若x是Rⁿ 中向量,则A与x的积,记为Ax,就是A的各列以x中对应元素为权的线性组合，即 [ Ax=egin{bmatrix} a_1 & a_2 & dots & a_n end{bmatrix} egin{bmatrix} x_1 \ x_2 \ vdots \ x_n end{bmatrix} = x_1a_1 + x_2a_2 + dots + x_n a_n ]

称Ax=b这样的方程为矩阵方程

至此，线性方程组有三种不同但彼此等价的观点:作为矩阵方程、作为向量方程或作为线性方程组。都用行化简算法来化简增广矩阵来解。

主对角线上元素为1，其他位置元素为0，这个矩阵称为单位矩阵。记为I

1.4 线性方程组的解集

如果线性方程组可以写成Ax=0的形式，则称为齐次的，A是mxn矩阵，0是Rⁿ 中的零向量。这样的方程组至少有一个解，即x=0,这个解称为它的平凡解。

齐次方程Ax=0有非平凡解，当且仅当方程至少有一个自由变量。

测试

线性代数学习笔记二

线性代数学习笔记二

目录

1 线性方程组

1.1 行化简与阶梯型矩阵

1.1.1 主元位置

1.1.2 线性方程组的解

1.2 向量方程

1.3 矩阵方程

1.4 线性方程组的解集