Linear Algebra Study Note

Views: 122
Wrote on May 25, 2020, 10:06 p.m.

Customized Linear Algebra Calculation Module

# _global.py
# -*- coding: utf-8 -*-
"""
Created by @南卡统爷 on 5/19/20 1:41 PM.
"""
__author__ = '@JiayuWang'

EPSILON = 1e-8


def is_zero(x):
    return abs(x) < EPSILON

def is_equal(a, b):
    return abs(a - b) < EPSILON
# Vector.py
# -*- coding: utf-8 -*-
"""
Created by @南卡统爷 on 5/18/20 9:55 PM.
"""
__author__ = '@JiayuWang'

import math

from ._global import is_zero, is_equal


class Vector:

    def __init__(self, lst):
        self._values = list(lst)

    @classmethod
    def zero(cls, dim):
        """返回一个dim维的零向量"""
        return cls([0] * dim)

    def norm(self):
        """返回向量的模(长度)"""
        return math.sqrt(sum(e ** 2 for e in self))

    def normalize(self):
        """返回向量的单位向量"""
        if is_zero(self.norm()):
            raise ZeroDivisionError("Normalize error! Norm is zero.")
        return 1 / self.norm() * Vector(self)  # 也可以写成 Vect(self._values) 但因为定义了 __iter__ 所以可以直接用 self
        # return Vector(self._values) / self.norm()

    def underlying_list(self):
        """返回向量的底层列表"""
        return self._values[:]

    def __add__(self, another):
        """向量加法,返回结果向量"""
        assert len(self) == len(another), \
            "Error in adding. Length of vectors must be equal"
        return Vector([a + b for a, b in zip(self, another)])

    def __sub__(self, another):
        """向量减法,返回结果向量"""
        assert len(self) == len(another), \
            "Error in adding. Length of vectors must be equal"
        return Vector([a - b for a, b in zip(self, another)])

    def dot(self, another):
        """向量点乘,返回结果标量"""
        assert len(self) == len(another), \
            "Error in adding. Length of vectors must be equal"
        return sum(a * b for a, b in zip(self, another))

    def __mul__(self, k):
        """向量乘法,返回结果向量: self * k """
        return Vector([k * e for e in self])

    def __rmul__(self, k):
        """向量乘法,返回结果向量: k * self  """
        return self * k

    def __truediv__(self, k):
        """返回数量除法的结果向量: self / k """
        return (1 / k) * self

    def __pos__(self):
        """返回向量取正的结果向量"""
        return 1 * self

    def __neg__(self):
        """返回向量取负的结果向量"""
        return -1 * self

    def __eq__(self, other):
        """返回向量是否相等"""
        other_list = other.underlying_list()
        if(len(other_list) != len(self._values)):
            return False
        return all(is_equal(x, y) for x, y in zip(self._values, other_list))

    def __neq__(self, other):
        """返回向量是否不等"""
        return not(self == other)

    def __iter__(self):
        """返回向量的迭代器"""
        return self._values.__iter__()

    def __getitem__(self, index):
        """取向量的第index个元素"""
        return self._values[index]

    def __len__(self):
        """返回向量长度(有多少个元素)"""
        return len(self._values)

    def __repr__(self):
        return "Vector({})".format(self._values)

    def __str__(self):
        return "({})".format(", ".join(str(e) for e in self._values))
# Matrix.py
# -*- coding: utf-8 -*-

__author__ = '@JiayuWang'

from .Vector import Vector


class Matrix:

    def __init__(self, list2d):

        if isinstance(list2d[0], list):
            self._values = [row[:] for row in list2d]
        elif isinstance(list2d[0], Vector):
            self._values = [row.underlying_list() for row in list2d]

    @classmethod
    def zero(cls, r, c):
        """返回一个r行c列的零矩阵"""
        return cls([[0]*c for _ in range(r)])

    @classmethod
    def identity(cls, n):
        """返回一个n行n列的零矩阵"""
        m = [[0]*n for _ in range(n)]
        for i in range(n):
            m[i][i] = 1;
        return cls(m)


    def T(self):
        """返回矩阵的转置矩阵"""
        return Matrix([[e for e in self.col_vector(i)] for i in range(self.col_num())])

    def __add__(self, another):
        """返回两个矩阵的相加结果"""
        assert self.shape() == another.shape(), \
            "Error in adding. Shape of matrix must be the same."
        return Matrix([[a+b for a,b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())])

    def __sub__(self, another):
        """返回两个矩阵的相减结果"""
        assert self.shape() == another.shape(), \
            "Error in adding. Shape of matrix must be the same."
        return Matrix([[a-b for a,b in zip(self.row_vector(i), another.row_vector(i))] for i in range(self.row_num())])

    def __mul__(self, k):
        """返回矩阵的数量乘法: self * k"""
        return Matrix([[e*k for e in self.row_vector(i)] for i in range(self.row_num())])

    def __rmul__(self, k):
        """返回矩阵的数量乘法: k * self"""
        return self * k

    def dot(self, another):
        """返回矩阵乘法的结果"""
        if isinstance(another, Vector):
            # 矩阵和向量的乘法
            assert self.col_num() == len(another), \
                "Error in Matrix-Vector Multiplication"
            return Vector([self.row_vector(i).dot(another) for i in range(self.row_num())])

        if isinstance(another, Matrix):
            # 矩阵和矩阵的乘法
            assert self.col_num() == another.row_num(), \
                "Error in Matrix-Matrix Multiplication"
            return Matrix([[self.row_vector(i).dot(another.col_vector(j)) for j in range(another.col_num())] for i in range(self.row_num())])

    def __truediv__(self, k):
        """返回矩阵的数量除法: self / k"""
        return (1/k) * self

    def __pos__(self):
        """返回矩阵的取正结果"""
        return 1*self

    def __neg__(self):
        """返回矩阵的取负结果"""
        return -1*self

    def row_vector(self, index):
        """返回矩阵的第index个行向量"""
        return Vector(self._values[index])

    def col_vector(self, index):
        """返回矩阵的第index个列向量"""
        return Vector([row[index] for row in self._values])
    def __getitem__(self, pos):
        """返回矩阵pos位置的元素"""
        r,c=pos
        return self._values[r][c]

    def size(self):
        """返回矩阵的元素个数"""
        r, c = self.shape()
        return r*c

    def shape(self):
        """返回矩阵的形状:(行数,列数)"""
        return len(self._values), len(self._values[0])

    def row_num(self):
        """返回矩阵的行数"""
        return self.shape()[0]

    __len__ = row_num

    def col_num(self):
        """返回矩阵的列数"""
        return self.shape()[1]

    def __repr__(self):
        return "Matrix({})".format(self._values)

    __str__ = __repr__
# LU.py
# -*- coding: utf-8 -*-
"""
Created by @南卡统爷 on 5/24/20 3:52 PM.
"""
__author__ = '@JiayuWang'

from .Matrix import Matrix
from .Vector import Vector
from ._global import is_zero


def lu(matrix):

    assert matrix.row_num() == matrix.col_num(), "Matrix must be a square matrix"

    n = matrix.row_num()
    A = [matrix.row_vector(i) for i in range(n)]
    L = [[1.0 if i == j else 0.0 for i in range(n)] for j in range(n)]

    for i in range(n):

        # 看A[i][i]位置是否可以是主元
        if is_zero(A[i][i]):
            return None, None
        else:
            for j in range(i+1, n):
                p = A[j][i] / A[i][i]
                A[j] = A[j] - p * A[i]
                L[j][i] = p

    return Matrix(L), Matrix([A[i].underlying_list() for i in range(n)])
# LinearSystem.py
# -*- coding: utf-8 -*-
"""
Created by @南卡统爷 on 5/23/20 1:19 AM.
"""
__author__ = '@JiayuWang'

from .Matrix import Matrix
from .Vector import Vector
from ._global import is_zero


class LinearSystem:

    def __init__(self, A, b=None):
        """A是增广矩阵左侧的系数矩阵,b是右侧的结果的列向量"""
        assert b is None or A.row_num() == len(b), "row number of A must be equal to the length of b"
        self._m = A.row_num()
        self._n = A.col_num()

        if b is None:
            self.Ab = [A.row_vector(i) for i in range(self._m)]
        if isinstance(b, Vector):
            self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]]) for i in range(self._m)]
        if isinstance(b, Matrix):
            self.Ab = [Vector(A.row_vector(i).underlying_list() + b.row_vector(i).underlying_list()) for i in
                       range(self._m)]

        self.pivots = []

    def _max_row(self, index_i, index_j, n):
        """找出指定列的数字最大的那行的行号"""
        best, ret = self.Ab[index_i][index_j], index_i
        for i in range(index_i + 1, n):
            if self.Ab[i][index_j] > best:
                best, ret = self.Ab[i][index_j], i
        return ret

    def _forward(self):

        i, k = 0, 0
        while i < self._m and k < self._n:
            # 看Ab[i][i]位置是否可以是主元
            max_row = self._max_row(i, k, self._m)
            self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i]

            if is_zero(self.Ab[i][k]):
                k += 1
            else:
                # 将主元归为一
                self.Ab[i] = self.Ab[i] / self.Ab[i][k]
                for j in range(i + 1, self._m):
                    self.Ab[j] = self.Ab[j] - self.Ab[j][k] * self.Ab[i]
                self.pivots.append(k)
                i += 1

    def _backward(self):

        n = len(self.pivots)
        for i in range(n - 1, -1, -1):
            k = self.pivots[i]
            # Ab[i][k]为主元
            for j in range(i - 1, -1, -1):
                self.Ab[j] = self.Ab[j] - self.Ab[j][i] * self.Ab[i]

    def gauss_jordan_elimination(self):
        """如果有解,返回True;如果没有解,返回False"""

        self._forward()
        self._backward()

        for i in range(len(self.pivots), self._m):
            if not is_zero(self.Ab[i][-1]):
                return False
        return True

    def fancy_print(self):

        for i in range(self._m):
            print(" ".join(str(self.Ab[i][j]) for j in range(self._n)), end=" ")
            print("|", self.Ab[i][-1])


def inv(A):

    if A.row_num() != A.col_num():
        return None

    n = A.row_num()
    ls = LinearSystem(A, Matrix.identity(n))
    if not ls.gauss_jordan_elimination():
        return None

    invA = [[row[i] for i in range(n, 2*n)] for row in ls.Ab]
    return Matrix(invA)


def rank(A):
    """求矩阵的秩"""
    ls = LinearSystem(A)
    ls.gauss_jordan_elimination()
    zero = Vector.zero(A.col_num())
    return sum([row != zero for row in ls.Ab])
# GramSchmidtProcess.py
# -*- coding: utf-8 -*-
"""
Created by @南卡统爷 on 5/26/20 1:29 AM.
"""
__author__ = '@JiayuWang'

from .Vector import Vector
from .Matrix import Matrix
from .LinearSystem import rank

def gram_schmidt_process(basis):
    #把所有的向量按照行的方式排列在一起然后高斯消元,只要结果没有零行就说明线性无关
    matrix = Matrix(basis)
    assert rank(matrix) == len(basis)

    res = [basis[0]]
    for i in range(1, len(basis)):
        p = basis[i]
        for r in res:
            p = p - basis[i].dot(r) / r.dot(r) * r
        res.append(p)
    return res


def qr(A):

    assert A.row_num() == A.col_num(), "A must be square"

    basis = [A.col_vector(i) for i in range(A.col_num())]
    P = gram_schmidt_process(basis)

    # 基于标准正交基
    Q = Matrix([v/v.norm() for v in P]).T()
    R = Q.T().dot(A)

    return Q, R