GMRES算法原理

• 第一步，格拉姆-施密特(Gram-Schmidt)过程计算高维投影

具体请参考线性代数整理(二) ，这里不再赘述

• 第二步，Krylov子空间及Arnoldi算法

设v是非零n维向量，A是n阶方阵，向量组生成的向量空间称为k维Krylov子空间(包含在(n维欧几里得空间)中)，记作，即，生成Krylov子空间的向量组称为Krylov向量组。这里假定Krylov向量组是线性无关的。

给定Krylov向量组(k维)，求一组与之等价的规范正交向量组。使用格拉姆-施密特过程求解便可。Arnoldi在研究非对称矩阵的特征值问题时，利用Krylov向量组的特殊结构，给出了格拉姆-施密特算法的一种变体算法，现在称为Arnoldi算法。关于特征值的具体内容请参考线性代数整理(三)

假设k=4，此时Krylov向量组为:，首先将v单位化，，然后用v1替换原向量组的v，得

该向量组与原向量组等价。Av1在v1上的投影记为，此处表示v1•Av1/||v1||，因为v1是单位向量，只表示方向，所以为投影向量，则从Av1中减去它在v1上的投影向量，得，v‘2是v1的正交向量，单位化后得，v2是v1的标准正交向量。用v2替换向量组中的Av1,得

这里替换后，由v1、Av1的这对非正交基变成了v1、v2的一对标准正交基，它们表示的是同一个平面。后面的就是代入替换，得。

可以证明向量组与等价，即它们可以相互线性相关表示。即中的每个向量可由中若干向量的线性组合表示。

证明：由以及，可得，或者写作

即Av1可由中的前两个向量的线性组合表示。在式两端左乘A，得

结合可知，可由中前三个向量的线性组合表示(因为Av1可以分解到中)。在式两端左乘A，得

结合可知，可由中的前四个向量的线性组合表示。

以上证明可由线性表示。

反过来，由可得，或者写作

即v2可由中的第一、二向量的线性组合表示，在式的两端左乘A，得

即Av2可由中的第二、三向量的线性组合表示，在式两端左乘A，得

即可由中的第三、四向量的线性组合表示，这就证明了可由向量组线性表示。

现在由出发，计算v3。从Av2中减去它在v1和v2上的投影向量，得

，v‘3是v1、v2组成平面的正交向量。单位化后得，v3是v1、v2组成平面的标准正交向量，用v3替换中的Av2，得向量组

这里替换后，由v1、v2、Av2的这组三维空间的非正交基变成了v1、v2、v3的三维空间的一组标准正交基，它们表示的是同一个三维空间。

与是等价向量组。现在由出发，计算v4，从Av3中减去它在v1、v2和v3上的投影向量，得

，v'4是与v1、v2、v3所构成的三维空间的正交向量。单位化后得，v4是与v1、v2、v3所构成的三维空间的标准正交向量。至此完成了k=4时的Arnoldi算法步骤，得到了标准正交向量组v1、v2、v3、v4。

这里跟格拉姆-施密特过程的区别只是把任意的一组基变成了的一组基。

伪代码

v[1] = v / ||v||      //v单位化
for (int j = 1;j < k;j++) {    //k维子空间
    w[j] = Av[j]               //生成下一个向量
    for (int i = 1;i <= j;i++) {
        h[i][j] = (w[j],v[i])   //获取Av[j]在v[i]上的投影向量的模
    }
    for (int i = 1;i <= j;i++) {
        w[j] -= h[i][j]v[i]    //获取j+1维空间与j维体的正交向量w[j]
    }
    h[j+1][j] = ||w[j]||       //获取该正交向量w[j]的模
    if (h[j+1][j] == 0) {
        break
    }
    v[j+1] = w[j] / h[j+1][j]  //w[j]单位化
}

最后得到的v[1]、v[2].......v[k]就是的一组标准正交基。此处(Av[j],v[i])表示v[i]•Av[j]/||v[i]||，表示的就是Av[j]在v[i]投影向量的模。

在上面的代码中，我们可以看到

h[j+1][j] = ||w[j]||       //获取该正交向量w[j]的模

即将每个 w[j] (j=1.....k)的模记录为h[j+1][j]，所有算得的构成一个(k+1)*k的上海森堡(Hessenberg)矩阵，对k=4的情形有：

这让我们想起了矩阵的QR分解，有关QR分解的内容请参考线性代数整理(二) ，这里被分解的矩阵就是k维Krylov子空间的这组基所组成的矩阵，分解后得到Q标准正交矩阵即为v1,v2,v3,...vk所组成的矩阵，R是一个上三角矩阵，虽然海森堡矩阵并不是一个上三角矩阵，但是是一个类上三角矩阵，这是因为所组成的矩阵并不是一个方阵，因为k是小于n的。

观察Arnoldi算法中的式子，即伪代码中的

w[j] = Av[j]               //生成下一个向量
for (int i = 1;i <= j;i++) {
    w[j] -= h[i][j]v[i]    //获取j+1维空间与j维体的正交向量w[j]
}

 及，即伪代码中的

v[j+1] = w[j] / h[j+1][j]  //w[j]单位化

立刻得到关系式，或写成矩阵形式：

                    (1)

将分成左右两块，分成上下两块，其中和分别表示的前k行和最后一行(的前k-1个元素都为零)，则有，在该式两端左乘，并注意到正交性，立即有

                             (2)

(这里因为是标准正交矩阵，所以就是的逆。)

Java代码样例

更新一下向量类

@Getter
@AllArgsConstructor
public class Vector<T extends Number> implements Attribute<T>,Cloneable {
    private T[] values;

    @SuppressWarnings("unchecked")
    public static <T extends Number> Attribute<T> zero(int dim,Addr<T> addr) {
        T[] values = (T[])new Number[dim];
        Arrays.fill(values,addr.zero());
        return new Vector<>(values);
    }

    @Override
    public String toString() {
        return "Vector{" +
                "values=" + Arrays.toString(values) +
                '}';
    }

    @Override
    @SuppressWarnings("unchecked")
    public Vector clone() throws CloneNotSupportedException {
        T[] values = Arrays.copyOf(this.values,this.len());
        return new Vector(values);
    }

    @Override
    @SuppressWarnings("unchecked")
    public boolean equals(Attribute<T> another,Addr<T> addr) {
        if (this.len() != another.len()) {
            return false;
        }
        for (int i = 0; i < len(); i++) {
            if (!addr.equals(this.get(i),another.get(i))) {
                return false;
            }
        }
        return true;
    }

    @Override
    public int len() {
        return values.length;
    }

    @Override
    public T get(int index) {
        if (index >= len() || index < 0) {
            throw new IllegalArgumentException("索引超出范围");
        }
        return values[index];
    }

    @Override
    public void set(int index, T value) {
        if (index >= len() || index < 0) {
            throw new IllegalArgumentException("索引超出范围");
        }
        this.values[index] = value;
    }

    @Override
    @SuppressWarnings("unchecked")
    public Attribute<T> add(Attribute<T> another,Addr<T> addr) {
        if (this.len() != another.len()) {
            throw new IllegalArgumentException("加法错误，两个向量的维度必须相等");
        }
        T[] newValues = (T[])new Number[this.len()];
        for (int i = 0; i < this.len(); i++) {
            newValues[i] = addr.add(this.get(i),another.get(i));
        }
        return new Vector<>(newValues);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Attribute<T> sub(Attribute<T> another, Addr<T> addr) {
        if (this.len() != another.len()) {
            throw new IllegalArgumentException("减法错误，两个向量的维度必须相等");
        }
        T[] newValues = (T[])new Number[this.len()];
        for (int i = 0; i < this.len(); i++) {
            newValues[i] = addr.sub(this.get(i),another.get(i));
        }
        return new Vector<>(newValues);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Attribute<T> div(T k, Addr<T> addr) {
        T[] newValues = (T[])new Number[this.len()];
        for (int i = 0; i < this.len(); i++) {
            newValues[i] = addr.div(this.get(i),k);
        }
        return new Vector<>(newValues);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Attribute<T> mul(T k, Addr<T> addr) {
        T[] newValues = (T[])new Number[this.len()];
        for (int i = 0; i < this.len(); i++) {
            newValues[i] = addr.mul(k,this.get(i));
        }
        return new Vector<>(newValues);
    }

    @Override
    public double norm(Addr<T> addr) {
        T total = addr.zero();
        for (int i = 0; i < this.len(); i++) {
            total = addr.add(total,addr.square(this.get(i)));
        }
        return addr.prescription(total);
    }

    @Override
    public Attribute<Double> normalize(Addr<T> addr) {
        Double[] newValues = new Double[this.len()];
        double norm = this.norm(addr);
        for (int i = 0; i < this.len(); i++) {
            newValues[i] = addr.div(this.get(i),norm);
        }
        return new Vector<>(newValues);
    }

    @Override
    public T dot(Attribute<T> another, Addr<T> addr) {
        if (this.len() != another.len()) {
            throw new IllegalArgumentException("点乘错误，两个向量的维度必须相等");
        }
        T total = addr.zero();
        for (int i = 0; i < this.len(); i++) {
            total = addr.add(total,addr.mul(this.get(i),another.get(i)));
        }
        return total;
    }

    @Override
    public double degree(Attribute<T> another, Addr<T> addr) {
        T dot = this.dot(another,addr);
        double selfNorm = this.norm(addr);
        double anotherNorm = another.norm(addr);
        double productNorm = selfNorm * anotherNorm;
        return addr.toDegree(addr.acos(addr.div(dot,productNorm)));
    }

    @Override
    public Attribute<Double> project(Attribute<T> another, Addr<T> addr) {
        T dot = this.dot(another,addr);
        double selfNorm = this.norm(addr);
        Attribute<Double> unitVector = this.normalize(addr);
        double d = addr.div(dot, selfNorm);
        return mul(d,unitVector);
    }

    private Attribute<Double> mul(double d,Attribute<Double> v) {
        Double[] values = new Double[v.len()];
        for (int i = 0; i < v.len(); i++) {
            values[i] = d * v.get(i);
        }
        return new Vector<>(values);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Attribute<Double> rectangular(Attribute<T> another, Addr<T> addr) {
        Attribute<Double> att = (Attribute<Double>) another;
        Addr<Double> add = (Addr<Double>) addr;
        return att.sub(this.project(another,addr),add);
    }
}

@AllArgsConstructor
public class Matrix<T extends Number> implements Ranks<T>,Cloneable {
    private T[][] values;

    @Override
    public String toString() {
        StringBuilder builder = new StringBuilder();
        for (int i = 0; i < values.length - 1; i++) {
            builder.append(Arrays.toString(values[i]) + ",");
        }
        builder.append(Arrays.toString(values[values.length - 1]));
        return "Matrix{" +
                "values=" + "[" + builder.toString() + "]" +
                '}';
    }

    /**
     * 创建一个零矩阵
     * @param row
     * @param col
     * @param addr
     * @param <T>
     * @return
     */
    @SuppressWarnings("unchecked")
    public static <T extends Number> Ranks<T> zero(int row,int col,Addr<T> addr) {
        T[][] values = (T[][])new Number[row][col];
        for (int i = 0; i < row; i++) {
            Arrays.fill(values[i],addr.zero());
        }
        return new Matrix<>(values);
    }

    /**
     * 创建一个单位矩阵
     * @param n
     * @param <T>
     * @return
     */
    @SuppressWarnings("unchecked")
    public static <T extends Number> Ranks<T> identity(int n,Addr<T> addr) {
        T[][] values = (T[][])new Number[n][n];
        for (int i = 0; i < n; i++) {
            Arrays.fill(values[i],addr.zero());
        }
        for (int i = 0; i < n; i++) {
            values[i][i] = addr.one();
        }
        return new Matrix<>(values);
    }

    @Override
    public int[] shape() {
        int[] res = new int[2];
        res[0] = values.length;
        res[1] = values[0].length;
        return res;
    }

    @Override
    public int rowNum() {
        return shape()[0];
    }

    @Override
    public int colNum() {
        return shape()[1];
    }

    @Override
    public int size() {
        int row = rowNum();
        int col = colNum();
        return row * col;
    }

    @Override
    public T get(int pos) {
        if (pos >= size() || pos < 0)  {
            throw new IllegalArgumentException("位置超出索引范围");
        }
        int row = pos / colNum();
        int col = pos % colNum();
        return values[row][col];
    }

    @Override
    public T get(Tuple<Integer, Integer> pos) {
        if ((int) pos.getField(0) >= rowNum() || (int) pos.getField(0) < 0 ||
                (int) pos.getField(1) >= colNum() || (int) pos.getField(1) < 0) {
            throw new IllegalArgumentException("位置超出索引范围");
        }
        return values[(int) pos.getField(0)][(int) pos.getField(1)];
    }

    @Override
    public T get(int[] pos) {
        if (pos.length != 2) {
            throw new IllegalArgumentException("位置长度超出范围");
        }
        if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) {
            throw new IllegalArgumentException("位置超出索引范围");
        }
        return values[pos[0]][pos[1]];
    }

    @Override
    public void set(int[] pos, T value) {
        if (pos.length != 2) {
            throw new IllegalArgumentException("位置长度超出范围");
        }
        if (pos[0] >= rowNum() || pos[0] < 0 || pos[1] >= colNum() || pos[1] < 0) {
            throw new IllegalArgumentException("位置超出索引范围");
        }
        values[pos[0]][pos[1]] = value;
    }

    @Override
    public Attribute<T> rowVector(int index) {
        if (index >= rowNum() || index < 0) {
            throw new IllegalArgumentException("索引超出范围");
        }
        return new Vector<>(values[index]);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Attribute<T> colVector(int index) {
        if (index >= colNum() || index < 0) {
            throw new IllegalArgumentException("索引超出范围");
        }
        T[] values = (T[])new Number[rowNum()];
        for (int i = 0; i < rowNum(); i++) {
            values[i] = this.values[i][index];
        }
        return new Vector<>(values);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Ranks<T> add(Ranks<T> another,Addr<T> addr) {
        if (!Arrays.equals(this.shape(),another.shape())) {
            throw new IllegalArgumentException("加法错误，两个矩阵的形状必须相同");
        }
        T[][] values = (T[][])new Number[this.rowNum()][this.colNum()];
        for (int i = 0; i < rowNum(); i++) {
            for (int j = 0; j < colNum(); j++) {
                values[i][j] = addr.add(this.get(new int[]{i,j}),another.get(new int[]{i,j}));
            }
        }
        this.values = values;
        return this;
    }

    @Override
    @SuppressWarnings("unchecked")
    public Ranks<T> sub(Ranks<T> another, Addr<T> addr) {
        if (!Arrays.equals(this.shape(),another.shape())) {
            throw new IllegalArgumentException("减法错误，两个矩阵的形状必须相同");
        }
        T[][] values = (T[][])new Number[this.rowNum()][this.colNum()];
        for (int i = 0; i < rowNum(); i++) {
            for (int j = 0; j < colNum(); j++) {
                values[i][j] = addr.sub(this.get(new int[]{i,j}),another.get(new int[]{i,j}));
            }
        }
        this.values = values;
        return this;
    }

    @Override
    @SuppressWarnings("unchecked")
    public Ranks<T> mul(T k, Addr<T> addr) {
        T[][] values = (T[][])new Number[this.rowNum()][this.colNum()];
        for (int i = 0; i < rowNum(); i++) {
            for (int j = 0; j < colNum(); j++) {
                values[i][j] = addr.mul(k,this.get(new int[]{i,j}));
            }
        }
        this.values = values;
        return this;
    }

    @Override
    public Ranks<Double> div(double k, Addr<T> addr) {
        Double[][] values = new Double[this.rowNum()][this.colNum()];
        for (int i = 0; i < rowNum(); i++) {
            for (int j = 0; j < colNum(); j++) {
                values[i][j] = addr.div(this.get(new int[]{i,j}),k);
            }
        }
        return new Matrix<>(values);
    }

    @Override
    public Ranks<T> pos(Addr<T> addr,T one) {
        return mul(one,addr);
    }

    @Override
    public Ranks<T> neg(Addr<T> addr,T negOne) {
        return mul(negOne,addr);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Attribute<T> dot(Attribute<T> another, Addr<T> addr) {
        if (this.colNum() != another.len()) {
            throw new IllegalArgumentException("矩阵列数需要与向量的长度相等");
        }
        T[] values = (T[])new Number[this.rowNum()];
        for (int i = 0; i < this.rowNum(); i++) {
            values[i] = this.rowVector(i).dot(another,addr);
        }
        return new Vector<>(values);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Ranks<T> dot(Ranks<T> another, Addr<T> addr) {
        if (this.colNum() != another.rowNum()) {
            throw new IllegalArgumentException("当前矩阵列数必须等于点乘矩阵行数");
        }
        T[][] values = (T[][]) new Number[this.rowNum()][another.colNum()];
        for (int i = 0; i < this.rowNum(); i++) {
            for (int j = 0; j < another.colNum(); j++) {
                values[i][j] = this.rowVector(i).dot(another.colVector(j),addr);
            }
        }
        return new Matrix<>(values);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Ranks<T> trans() {
        T[][] values = (T[][])new Number[this.colNum()][this.rowNum()];
        for (int i = 0; i < this.colNum(); i++) {
            for (int j = 0; j < this.rowNum(); j++) {
                values[i][j] = this.values[j][i];
            }
        }
        return new Matrix<>(values);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Ranks<T> invFunc(Addr<T> addr) {
        if (this.rowNum() != this.colNum()) {
            return null;
        }
        int n = this.rowNum();
        FuncGroup<T> ls = new LinearSystemFunc<>(this,Matrix.identity(n,addr));
        if (!ls.gaussJordanElimination(addr)) {
            return null;
        }
        T[][] values = (T[][])new Number[this.rowNum()][((LinearSystemFunc)ls).getAb().colNum() - this.colNum()];
        for (int i = 0; i < this.rowNum(); i++) {
            for (int j = 0; j < ((LinearSystemFunc)ls).getAb().colNum() - this.colNum(); j++) {
                values[i][j] = (T)((LinearSystemFunc)ls).getAb().get(new int[]{i,((LinearSystemFunc)ls).getAb().colNum() - this.colNum() + j});
            }
        }
        return new Matrix<>(values);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Ranks<T> invMatrix(Addr<T> addr) {
        if (this.rowNum() != this.colNum()) {
            return null;
        }
        int n = this.rowNum();
        FuncGroup<T> ls = new LinearSystemMatrix<>(this,Matrix.identity(n,addr));
        if (!ls.gaussJordanElimination(addr)) {
            return null;
        }
        T[][] values = (T[][])new Number[this.rowNum()][((LinearSystemMatrix)ls).getAb().colNum() - this.colNum()];
        for (int i = 0; i < this.rowNum(); i++) {
            for (int j = 0; j < ((LinearSystemMatrix)ls).getAb().colNum() - this.colNum(); j++) {
                values[i][j] = (T)((LinearSystemMatrix)ls).getAb().get(new int[]{i,((LinearSystemMatrix)ls).getAb().colNum() - this.colNum() + j});
            }
        }
        return new Matrix<>(values);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Tuple<Ranks<T>,Ranks<T>> lu(Addr<T> addr) {
//        if (this.rowNum() != this.colNum()) {
//            throw new IllegalArgumentException("矩阵必须为方阵");
//        }
        int n = this.rowNum();
        int l = this.colNum();
        if (l > n) {
            T[][] A = (T[][])new Number[n][l];
            //此处不能直接把this.values赋给A，否则会改变原矩阵的values
            for (int i = 0; i < n; i++) {
                A[i] = Arrays.copyOf(this.values[i],n);
            }
        }
        T[][] A = (T[][])new Number[n][n];
        //此处不能直接把this.values赋给A，否则会改变原矩阵的values
        for (int i = 0; i < n; i++) {
            A[i] = Arrays.copyOf(this.values[i],n);
        }
        //构建一个下三角矩阵的二维数组，先将其初始为一个单位矩阵
        T[][] L = (T[][])new Number[n][n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                if (i == j) {
                    L[i][j] = addr.one();
                }else {
                    L[i][j] = addr.zero();
                }
            }
        }
        for (int i = 0; i < n; i++) {
            //主元不能为0
            if (addr.equals(A[i][i],addr.zero())) {
                return null;
            }else {
                //将主元下面的元素变为0
                for (int j = i + 1; j < n; j++) {
                    //获取初等矩阵逆矩阵的系数
                    T p = addr.div(A[j][i],A[i][i]);
                    for (int k = 0; k < n; k++) {
                        T u = A[i][k];
                        A[j][k] = addr.sub(A[j][k],addr.mul(p,u));
                    }
                    //将初等矩阵逆矩阵的系数填充下三角矩阵
                    L[j][i] = p;
                }
            }
        }
        return new Tuple<>(new Matrix(L),new Matrix(A));
    }

    @Override
    @SuppressWarnings("unchecked")
    public int rank(Addr<T> addr) {
        Attribute<T> zero = Vector.zero(this.rowNum(),addr);
        FuncGroup<T> ls = new LinearSystemFunc<>(this,zero);
        ls.gaussJordanElimination(addr);
        zero = Vector.zero(this.colNum() + 1,addr);
        int sum = 0;
        for (int i = 0; i < ((LinearSystem) ls).getAb().rowNum(); i++) {
            Attribute<T> row = ((LinearSystem) ls).getAb().rowVector(i);
            //找出不为零向量的行数，即矩阵的秩，空间维度
            if (!zero.equals(row,addr)) {
                sum++;
            }
        }
        return sum;
    }

    @Override
    @SuppressWarnings("unchecked")
    public List<Attribute<T>> rowBase(Addr<T> addr) {
        Attribute<T> zero = Vector.zero(this.rowNum(),addr);
        FuncGroup<T> ls = new LinearSystemFunc<>(this,zero);
        ls.gaussJordanElimination(addr);
        zero = Vector.zero(this.colNum() + 1,addr);
        List<Attribute<T>> rowBaseList = new ArrayList<>();
        for (int i = 0; i < ((LinearSystem) ls).getAb().rowNum(); i++) {
            Attribute<T> row = ((LinearSystem) ls).getAb().rowVector(i);
            //找出行最简形式不为零向量的向量
            if (!zero.equals(row,addr)) {
                T[] values = Arrays.copyOf((T[]) ((Vector) row).getValues(),this.colNum());
                Attribute<T> newRow = new Vector<>(values);
                rowBaseList.add(newRow);
            }
        }
        return rowBaseList;
    }

    @Override
    @SuppressWarnings("unchecked")
    public List<Attribute<T>> colBase(Addr<T> addr) {
        Attribute<T> zero = Vector.zero(this.rowNum(),addr);
        FuncGroup<T> ls = new LinearSystemFunc<>(this,zero);
        ls.gaussJordanElimination(addr);
        List<Integer> pivots = ((LinearSystem) ls).getPivots();
        List<Attribute<T>> colBaseList = new ArrayList<>();
        pivots.forEach(i -> colBaseList.add(this.colVector(i)));
        return colBaseList;
    }

    @Override
    public int dimensionNullSpace(Addr<T> addr) {
        return this.colNum() - rank(addr);
    }

    @Override
    @SuppressWarnings("unchecked")
    public List<Attribute<T>> nullBase(Addr<T> addr) {
        Attribute<T> zero = Vector.zero(this.rowNum(),addr);
        FuncGroup<T> ls = new LinearSystemFunc<>(this,zero);
        ls.gaussJordanElimination(addr);
        List<Integer> pivots = ((LinearSystem) ls).getPivots();
        List<Attribute<T>> zeroBaseList = new ArrayList<>();
        for (int i = 0; i < this.colNum(); i++) {
            //如果是自由列
            if (!pivots.contains(i)) {
                Attribute<T> temp = Vector.zero(this.colNum(),addr);
                Attribute<T> col = ((LinearSystem) ls).getAb().colVector(i);
                for (int j = 0; j < temp.len(); j++) {
                    if (j < col.len() && !addr.equals(col.get(j),addr.zero())) {
                        temp.set(j, addr.neg(col.get(j)));
                    }
                    if (i == j) {
                        temp.set(j,addr.one());
                    }
                }
                zeroBaseList.add(temp);
            }
        }
        return zeroBaseList;
    }

    @Override
    public int dimensionLeftNullSpace(Addr<T> addr) {
        return this.rowNum() - rank(addr);
    }

    @Override
    @SuppressWarnings("unchecked")
    public List<Attribute<T>> leftNullBase(Addr<T> addr) {
        return this.trans().nullBase(addr);
    }

    /**
     * 根据空间的一组基生成空间的一组正交基
     * @param basis 空间的一组基
     * @param addr
     * @return
     */
    @SuppressWarnings("unchecked")
    public static <T extends Number> List<Attribute<Double>> gramSchmidtProcess(List<Attribute<T>> basis, Addr<T> addr) {
        T[][] values = (T[][]) new Number[basis.size()][basis.get(0).len()];
        for (int i = 0; i < basis.size(); i++) {
            for (int j = 0; j < basis.get(0).len(); j++) {
                values[i][j] = basis.get(i).get(j);
            }
        }
        Ranks<T> matrix = new Matrix<>(values);
        if (matrix.rank(addr) != basis.size()) {
            throw new IllegalArgumentException("列表中的向量并非空间的一组基");
        }
        List<Attribute<Double>> res = new ArrayList<>();
        res.add((Attribute<Double>) basis.get(0));
        Addr<Double> add = (Addr<Double>) addr;
        //获取每i个正交基的生成空间的正交向量
        //第一次获取任意基的第一个向量和第二个向量所生成的二维平面的第一个向量的正交向量
        //第二次获取任意基的第一个向量和获取的正交向量生成的二维空间的三维正交向量
        //第三次获取三个正交向量生成的三维空间的四维正交向量，以此类推
        //所有的正交向量最终构成正交基
        for (int i = 1; i < basis.size(); i++) {
            Attribute<Double> p = (Attribute<Double>) basis.get(i);
            for (Attribute<Double> r : res) {
                p = r.rectangular(p,add);
            }
            res.add(p);
        }
        return res;
    }

    /**
     * Arnoldi算法
     * @param A 初始非对称n阶方阵
     * @param v 初始非零n维向量
     * @param k Krylov子空间维度
     * @param addr
     * @return
     */
    @SuppressWarnings("unchecked")
    public static Tuple<List<Attribute<Double>>,Double[][]> arnoldi(Ranks<Double> A,Attribute<Double> v,int k,Addr<Double> addr) {
        //k维Krylov子空间的标准正交向量集
        List<Attribute<Double>> u = new ArrayList<>();
        //海森堡矩阵
        Double[][] h = new Double[k + 1][k];
        for (int i = 0; i < k + 1; i++) {
            Arrays.fill(h[i],addr.zero());
        }
        //正交向量集
        Attribute<Double>[] w = new Attribute[k];
        //将初始向量单位化放入标准正交向量集中
        u.add(v.normalize(addr));
        for (int j = 0; j < k; j++) {
            //生成Krylov子空间的基的下一个向量
            w[j] = A.dot(u.get(j),addr);
            for (int i = 0; i <= j; i++) {
                //获取Av[j]在v[i]上的投影向量的模，并构建海森堡矩阵
                h[i][j] = u.get(i).dot(A.dot(u.get(j), addr), addr) / u.get(i).norm(addr);
            }
            //获取j+1维空间与j维体的正交向量
            for (int i = 0; i <= j; i++) {
                w[j] = w[j].sub(u.get(i).mul(h[i][j],addr),addr);
            }
            //获取该正交向量的模
            h[j + 1][j] = w[j].norm(addr);
            if (addr.equals(h[j + 1][j],addr.zero())) {
                break;
            }
            //该正交向量单位化并放入标准正交向量集中
            u.add(w[j].normalize(addr));
        }
        return new Tuple<>(u,h);
    }

    /**
     * 根据空间的一组基生成空间的一组标准正交基
     * @param basis 空间的一组基
     * @param addr
     * @param <T>
     * @return
     */
    @SuppressWarnings("unchecked")
    public static <T extends Number> List<Attribute<Double>> orthonormalBasis(List<Attribute<T>> basis, Addr<T> addr) {
        List<Attribute<Double>> bas = gramSchmidtProcess(basis,addr);
        Addr<Double> add = (Addr<Double>) addr;
        List<Attribute<Double>> res = bas.stream().map(vec -> vec.normalize(add))
                .collect(Collectors.toList());
        return res;
    }

    /**
     * 由一组标准正交基生成一个标准正交矩阵
     * @param basis 标准正交基
     * @return
     */
    public static Ranks<Double> orthonormalMatrix(List<Attribute<Double>> basis) {
        Double[][] values = new Double[basis.get(0).len()][basis.size()];
        for (int i = 0; i < basis.get(0).len(); i++) {
            for (int j = 0; j < basis.size(); j++) {
                values[i][j] = basis.get(j).get(i);
            }
        }
        return new Matrix<>(values);
    }

    @Override
    @SuppressWarnings("unchecked")
    public Tuple<Ranks<Double>,Ranks<Double>> qr(Addr<T> addr) {
        if (this.rowNum() != this.colNum()) {
            throw new RuntimeException("矩阵必须为一个方阵");
        }
        if (this.rank(addr) != this.colNum()) {
            throw new RuntimeException("矩阵的列向量必须线性无关");
        }
        List<Attribute<T>> basis = new ArrayList<>();
        for (int i = 0; i < this.colNum(); i++) {
            basis.add(this.colVector(i));
        }
        List<Attribute<Double>> sp = Matrix.orthonormalBasis(basis,addr);
        Ranks<Double> Q = Matrix.orthonormalMatrix(sp);
        Ranks<Double> source = (Ranks<Double>) this;
        Addr<Double> add = (Addr<Double>) addr;
        Ranks<Double> R = Q.trans().dot(source,add);
        return new Tuple<>(Q,R);
    }

    /**
     * 其他坐标系坐标转换标准坐标系坐标
     * @param otherBasis 其他坐标系的一组基
     * @param otherVector 其他坐标系的一个向量坐标
     * @param <T>
     * @return
     */
    public static <T extends Number> Attribute<T> otherToStandard(List<Attribute<T>> otherBasis,
                                                                  Attribute<T> otherVector,
                                                                  Addr<T> addr) {
        Ranks<T> matrix = buildColMatrix(otherBasis,addr);
        return matrix.dot(otherVector,addr);
    }

    /**
     * 标准坐标系坐标转换其他坐标系坐标
     * @param otherBasis 其他坐标系的一组基
     * @param vector 标准坐标系的一个向量坐标
     * @param addr
     * @param <T>
     * @return
     */
    public static <T extends Number> Attribute<T> standardToOther(List<Attribute<T>> otherBasis,
                                                                  Attribute<T> vector,
                                                                  Addr<T> addr) {
        Ranks<T> matrix = buildColMatrix(otherBasis,addr);
        return matrix.invFunc(addr).dot(vector,addr);
    }

    /**
     * 一个坐标系坐标转换另一个坐标系坐标
     * @param basis 第一个坐标系的一组基
     * @param anotherBasis 另一个坐标系的一组基
     * @param vector 第一个坐标系的一个向量坐标
     * @param addr
     * @param <T>
     * @return
     */
    public static <T extends Number> Attribute<T> oneToAnother(List<Attribute<T>> basis,
                                                               List<Attribute<T>> anotherBasis,
                                                               Attribute<T> vector,
                                                               Addr<T> addr) {
        Attribute<T> standardVector = otherToStandard(basis,vector,addr);
        return standardToOther(anotherBasis,standardVector,addr);
    }

    /**
     * 根据空间的一组基建立一个列向量矩阵
     * @param basis 空间的一组基
     * @param addr
     * @param <T>
     * @return
     */
    @SuppressWarnings("unchecked")
    private static <T extends Number> Ranks<T> buildColMatrix(List<Attribute<T>> basis,Addr<T> addr) {
        T[][] values = (T[][]) new Number[basis.get(0).len()][basis.size()];
        for (int i = 0; i < basis.get(0).len(); i++) {
            for (int j = 0; j < basis.size(); j++) {
                values[i][j] = basis.get(j).get(i);
            }
        }
        Ranks<T> matrix = new Matrix<>(values);
        if (matrix.rowNum() != matrix.colNum()) {
            throw new IllegalArgumentException("矩阵必须为一个方阵");
        }
        if (matrix.rank(addr) != basis.size()) {
            throw new IllegalArgumentException("列表中的向量并非空间的一组基");
        }
        return matrix;
    }

    @Override
    @SuppressWarnings("unchecked")
    public List<T> eigenvalue(Addr<T> addr) throws CloneNotSupportedException {
        if (this.rowNum() != this.colNum()) {
            throw new RuntimeException("当前矩阵必须为方阵");
        }
        List<T> res = new ArrayList<>();
        for (Double x = 0.0; x < 10.0; x++) {
            Ranks<T> I = Matrix.identity(this.colNum(),addr);
            Ranks<T> temp = this.clone();
            temp.sub(I.mul((T) x, addr), addr);
            Formula<T> determinant = new Determinant<>(temp);
            if (addr.equals(determinant.result(addr),addr.zero())) {
                res.add((T) x);
            }
        }
        return res;
    }

    @Override
    @SuppressWarnings("unchecked")
    public Matrix clone() throws CloneNotSupportedException {
        T[][] values = (T[][])new Number[this.rowNum()][this.colNum()];
        for (int i = 0; i < rowNum(); i++) {
            values[i] = Arrays.copyOf(this.values[i],this.colNum());
        }
        return new Matrix(values);
    }

    @SuppressWarnings("unchecked")
    public static void main(String[] args) throws CloneNotSupportedException {
        Ranks<Integer> matrix = new Matrix<>(new Integer[][]{{1,2,3},{3,4,5}});
        System.out.println(matrix);
        System.out.println(Arrays.toString(matrix.shape()));
        System.out.println(matrix.rowNum());
        System.out.println(matrix.colNum());
        System.out.println(matrix.size());
        System.out.println(matrix.get(4));
        System.out.println(matrix.get(new int[]{0,0}));
        System.out.println(matrix.rowVector(1));
        System.out.println(matrix.colVector(2));
        Ranks<Integer> matrix1 = new Matrix<>(new Integer[][]{{7,8,9},{10,20,30}});
        Addr<Integer> addr = new Addr<Integer>() {
            @Override
            public Integer zero() {
                return 0;
            }

            @Override
            public Integer one() {
                return 1;
            }

            @Override
            public Integer add(Integer lhs, Integer rhs) {
                return lhs + rhs;
            }

            @Override
            public Integer sub(Integer lhs, Integer rhs) {
                return lhs - rhs;
            }

            @Override
            public Integer mul(Integer k, Integer hs) {
                return k * hs;
            }

            @Override
            public Integer square(Integer hs) {
                return hs * hs;
            }

            @Override
            public double prescription(Integer hs) {
                return Math.sqrt(hs);
            }

            @Override
            public double div(Integer hs, double nhs) {
                return hs / nhs;
            }

            @Override
            public Integer div(Integer hs, Integer nhs) {
                return hs / nhs;
            }

            @Override
            public double acos(double hs) {
                return Math.acos(hs);
            }

            @Override
            public double toDegree(double hs) {
                return Math.toDegrees(hs);
            }

            @Override
            public int compair(Integer lhs, Integer rhs) {
                if (lhs - rhs > 0) {
                    return 1;
                }else if (lhs - rhs < 0) {
                    return -1;
                }else {
                    return 0;
                }
            }

            @Override
            public boolean equals(Integer lhs, Integer rhs) {
                return lhs == rhs;
            }

            @Override
            public Integer neg(Integer hs) {
                return -hs;
            }
        };
        System.out.println(matrix.trans());
        Ranks<Integer> copy = ((Matrix)matrix).clone();
        Attribute<Integer> vector = new Vector<>(new Integer[]{9,33,97});
        Ranks<Integer> matrix2 = new Matrix<>(new Integer[][]{{7,9},{11,33},{42,97}});
        //[(1,2,3),(3,4,5)]*(23,65,79)
        System.out.println(matrix.dot(vector,addr));
//        [(1,2,3),(3,4,5)]*[(7,9),(11,33),(42,97)]
        System.out.println(matrix.dot(matrix2,addr));
        Ranks<Integer> copy1 = ((Matrix)matrix).clone();
        System.out.println(matrix.add(matrix1,addr));
        System.out.println(matrix1.sub(copy,addr));
        System.out.println(copy.div(2,addr));
        System.out.println(copy.pos(addr,1));
        System.out.println(copy.neg(addr,-1));
        System.out.println(copy1.mul(3,addr));
        System.out.println(Matrix.zero(2,3,addr));
        Ranks<Integer> I = Matrix.identity(4, addr);
        System.out.println(I);
        Matrix<Integer> idenTest = new Matrix<>(new Integer[][]{{1,2,3,4},{5,6,7,8},{9,10,11,12},{13,14,15,16}});
        System.out.println(idenTest.dot(I,addr));
        System.out.println(I.dot(idenTest,addr));
        Addr<Double> addr1 = new Addr<Double>() {
            @Override
            public Double zero() {
                return 0.0;
            }

            @Override
            public Double one() {
                return 1.0;
            }

            @Override
            public Double add(Double lhs, Double rhs) {
                return lhs + rhs;
            }

            @Override
            public Double sub(Double lhs, Double rhs) {
                return lhs - rhs;
            }

            @Override
            public Double mul(Double k, Double hs) {
                return k * hs;
            }

            @Override
            public Double square(Double hs) {
                return hs * hs;
            }

            @Override
            public double prescription(Double hs) {
                return Math.sqrt(hs);
            }

            @Override
            public double div(Double hs, double nhs) {
                return hs / nhs;
            }

            @Override
            public Double div(Double hs, Double nhs) {
                return hs / nhs;
            }

            @Override
            public double acos(double hs) {
                return Math.acos(hs);
            }
            @Override
            public double toDegree(double hs) {
                return Math.toDegrees(hs);
            }

            @Override
            public int compair(Double lhs, Double rhs) {
                if (lhs - rhs > 0) {
                    return 1;
                }else if (lhs - rhs < 0) {
                    return -1;
                }else {
                    return 0;
                }
            }

            @Override
            public boolean equals(Double lhs, Double rhs) {
                double dis = 1e-8;
                if (Math.abs(lhs - rhs) < dis) {
                    return true;
                }
                return false;
            }

            @Override
            public Double neg(Double hs) {
                return -hs;
            }
        };
        Ranks<Double> A = new Matrix<>(new Double[][]{{1.0,2.0},{3.0,4.0}});
        System.out.println(A.invFunc(addr1));
        System.out.println(A.invMatrix(addr1));
        Ranks<Integer> B = new Matrix<>(new Integer[][]{{1,2,3},{4,5,6},{3,-3,5}});
        Tuple<Ranks<Integer>,Ranks<Integer>> lu = B.lu(addr);
        Ranks<Integer> L = lu.getField(0);
        System.out.println(L);
        Ranks<Integer> U = lu.getField(1);
        System.out.println(U);
        System.out.println(L.dot(U,addr));
        Ranks<Integer> N = new Matrix<>(new Integer[][]{{2,0,0},{-1,0,0},{0,0,1}});
        System.out.println(N.rank(addr));
        System.out.println(N.rowBase(addr));
        System.out.println(N.colBase(addr));
        System.out.println(N.dimensionNullSpace(addr));
        System.out.println(N.nullBase(addr));
        Ranks<Integer> M = new Matrix<>(new Integer[][]{{-1,2,3},{1,-4,-13},{-3,5,4}});
        System.out.println(M.rank(addr));
        System.out.println(M.rowBase(addr));
        System.out.println(M.colBase(addr));
        System.out.println(M.dimensionNullSpace(addr));
        System.out.println(M.nullBase(addr));
        Ranks<Double> F = new Matrix<>(new Double[][]{{1.0,2.0,3.0,4.0,5.0,6.0},
                {27.0,28.0,29.0,30.0,31.0,32.0},{15.0,16.0,17.0,18.0,19.0,20.0},
                {31.0,32.0,33.0,34.0,35.0,36.0}, {45.0,46.0,47.0,48.0,49.0,50.0}});
        System.out.println(F.rank(addr1));
        System.out.println(F.trans());
        System.out.println(F.trans().rank(addr1));
        System.out.println(F.rowBase(addr1));
        System.out.println(F.colBase(addr1));
        System.out.println(F.dimensionNullSpace(addr1));
        System.out.println(F.nullBase(addr1));
        System.out.println(F.dimensionLeftNullSpace(addr1));
        System.out.println(F.leftNullBase(addr1));
        List<Attribute<Double>> basis = Arrays.asList(new Vector<>(new Double[]{2.0,1.0}),
                new Vector<>(new Double[]{1.0,1.0}));
        System.out.println(Matrix.gramSchmidtProcess(basis,addr1));
        System.out.println(Matrix.orthonormalBasis(basis,addr1));
        Ranks<Double> orthonormalMatrix = Matrix.orthonormalMatrix(Matrix.orthonormalBasis(basis,addr1));
        System.out.println(orthonormalMatrix);
        System.out.println(orthonormalMatrix.dot(orthonormalMatrix.trans(),addr1));
        Ranks<Double> A1 = new Matrix<>(new Double[][]{{1.0,1.0,2.0},{1.0,1.0,0.0},
                {1.0,0.0,0.0}});
        Tuple<Ranks<Double>,Ranks<Double>> qr = A1.qr(addr1);
        Ranks<Double> Q = qr.getField(0);
        Ranks<Double> R = qr.getField(1);
        System.out.println(Q);
        System.out.println(R);
        System.out.println(Q.dot(R,addr1));
        List<Attribute<Double>> otherBasis = Arrays.asList(new Vector<>(new Double[]{4.0,1.0}),
                new Vector<>(new Double[]{2.0,3.0}));
        Attribute<Double> otherVector = new Vector<>(new Double[]{2.0,2.0});
        Attribute<Double> vectorStard = Matrix.otherToStandard(otherBasis,otherVector,addr1);
        System.out.println(vectorStard);
        System.out.println(Matrix.standardToOther(otherBasis,vectorStard,addr1));
        List<Attribute<Double>> anotherBasis = Arrays.asList(new Vector<>(new Double[]{1.0,3.0}),
                new Vector<>(new Double[]{2.0,4.0}));
        Attribute<Double> vectorAnother = Matrix.oneToAnother(otherBasis,anotherBasis,otherVector,addr1);
        System.out.println(vectorAnother);
        System.out.println(Matrix.otherToStandard(anotherBasis,vectorAnother,addr1));
        Ranks<Double> m = new Matrix<>(new Double[][]{{4.0,-2.0},{1.0,1.0}});
        System.out.println(m.eigenvalue(addr1));
        Ranks<Double> A2 = new Matrix<>(new Double[][]{{1.0,2.0},{3.0,4.0}});
        Attribute<Double> v = new Vector<>(new Double[]{2.0,2.0});
        Tuple<List<Attribute<Double>>,Double[][]> arnoldi = Matrix.arnoldi(A2, v, 2, addr1);
        List<Attribute<Double>> u = arnoldi.getField(0);
        Double[][] h = arnoldi.getField(1);
        System.out.println(u);
        for (int i = 0; i < h.length; i++) {
            System.out.println(Arrays.toString(h[i]));
        }
    }
}

运行结果

Matrix{values=[[1, 2, 3],[3, 4, 5]]}
[2, 3]
2
3
6
4
1
Vector{values=[3, 4, 5]}
Vector{values=[3, 5]}
Matrix{values=[[1, 3],[2, 4],[3, 5]]}
Vector{values=[366, 644]}
Matrix{values=[[155, 366],[275, 644]]}
Matrix{values=[[8, 10, 12],[13, 24, 35]]}
Matrix{values=[[6, 6, 6],[7, 16, 25]]}
Matrix{values=[[0.5, 1.0, 1.5],[1.5, 2.0, 2.5]]}
Matrix{values=[[1, 2, 3],[3, 4, 5]]}
Matrix{values=[[-1, -2, -3],[-3, -4, -5]]}
Matrix{values=[[3, 6, 9],[9, 12, 15]]}
Matrix{values=[[0, 0, 0],[0, 0, 0]]}
Matrix{values=[[1, 0, 0, 0],[0, 1, 0, 0],[0, 0, 1, 0],[0, 0, 0, 1]]}
Matrix{values=[[1, 2, 3, 4],[5, 6, 7, 8],[9, 10, 11, 12],[13, 14, 15, 16]]}
Matrix{values=[[1, 2, 3, 4],[5, 6, 7, 8],[9, 10, 11, 12],[13, 14, 15, 16]]}
Matrix{values=[[-1.9999999999999996, 0.9999999999999998],[1.4999999999999998, -0.49999999999999994]]}
Matrix{values=[[-1.9999999999999996, 0.9999999999999998],[1.4999999999999998, -0.4999999999999999]]}
Matrix{values=[[1, 0, 0],[4, 1, 0],[3, 3, 1]]}
Matrix{values=[[1, 2, 3],[0, -3, -6],[0, 0, 14]]}
Matrix{values=[[1, 2, 3],[4, 5, 6],[3, -3, 5]]}
2
[Vector{values=[1, 0, 0]}, Vector{values=[0, 0, 1]}]
[Vector{values=[2, -1, 0]}, Vector{values=[0, 0, 1]}]
1
[Vector{values=[0, 1, 0]}]
2
[Vector{values=[1, 0, 7]}, Vector{values=[0, 1, 5]}]
[Vector{values=[-1, 1, -3]}, Vector{values=[2, -4, 5]}]
1
[Vector{values=[-7, -5, 1]}]
2
Matrix{values=[[1.0, 27.0, 15.0, 31.0, 45.0],[2.0, 28.0, 16.0, 32.0, 46.0],[3.0, 29.0, 17.0, 33.0, 47.0],[4.0, 30.0, 18.0, 34.0, 48.0],[5.0, 31.0, 19.0, 35.0, 49.0],[6.0, 32.0, 20.0, 36.0, 50.0]]}
2
[Vector{values=[1.0, 0.0, -0.9999999999999993, -1.9999999999999998, -3.0, -3.9999999999999996]}, Vector{values=[0.0, 1.0, 1.9999999999999998, 3.0, 4.0, 5.0]}]
[Vector{values=[1.0, 27.0, 15.0, 31.0, 45.0]}, Vector{values=[2.0, 28.0, 16.0, 32.0, 46.0]}]
4
[Vector{values=[0.9999999999999993, -1.9999999999999998, 1.0, 0.0, 0.0, 0.0]}, Vector{values=[1.9999999999999998, -3.0, 0.0, 1.0, 0.0, 0.0]}, Vector{values=[3.0, -4.0, 0.0, 0.0, 1.0, 0.0]}, Vector{values=[3.9999999999999996, -5.0, 0.0, 0.0, 0.0, 1.0]}]
3
[Vector{values=[-0.4615384615384621, -0.5384615384615384, 1.0, 0.0, 0.0]}, Vector{values=[0.1538461538461533, -1.1538461538461537, 0.0, 1.0, 0.0]}, Vector{values=[0.6923076923076898, -1.692307692307692, 0.0, 0.0, 1.0]}]
[Vector{values=[2.0, 1.0]}, Vector{values=[-0.19999999999999996, 0.4]}]
[Vector{values=[0.8944271909999159, 0.4472135954999579]}, Vector{values=[-0.44721359549995787, 0.894427190999916]}]
Matrix{values=[[0.8944271909999159, -0.44721359549995787],[0.4472135954999579, 0.894427190999916]]}
Matrix{values=[[0.9999999999999999, 0.0],[0.0, 1.0000000000000002]]}
Matrix{values=[[0.5773502691896258, 0.4082482904638628, 0.7071067811865475],[0.5773502691896258, 0.4082482904638628, -0.7071067811865476],[0.5773502691896258, -0.8164965809277263, 0.0]]}
Matrix{values=[[1.7320508075688776, 1.1547005383792517, 1.1547005383792517],[-6.661338147750939E-16, 0.8164965809277256, 0.8164965809277256],[-1.1102230246251565E-16, -1.1102230246251565E-16, 1.414213562373095]]}
Matrix{values=[[0.9999999999999999, 0.9999999999999997, 1.9999999999999996],[1.0, 0.9999999999999999, -2.220446049250313E-16],[1.0000000000000007, 3.3306690738754696E-16, 3.3306690738754696E-16]]}
Vector{values=[12.0, 8.0]}
Vector{values=[1.9999999999999996, 2.0]}
Vector{values=[-15.999999999999995, 13.999999999999996]}
Vector{values=[11.999999999999998, 8.0]}
[2.0, 3.0]
[Vector{values=[0.7071067811865475, 0.7071067811865475]}, Vector{values=[-0.7071067811865476, 0.7071067811865474]}]
[5.0, 0.999999999999999]
[1.9999999999999996, -5.551115123125784E-16]
[0.0, 1.8054905945397061E-16]

Python代码

from .Vector import Vector
from ._global import is_zero

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 trans(self):
        # 返回矩阵的转置矩阵
        return Matrix([[e for e in self.col_vector(i)]
                       for i in range(self.col_num())])

    def inv_func(self):
        # 返回矩阵的逆
        from .LinearSystemFunc import LinearSystemFunc
        if self.row_num() != self.col_num():
            return None
        n = self.row_num()
        ls = LinearSystemFunc(self, 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 inv_matrix(self):
        # 返回矩阵的逆
        from .LinearSystemMatrix import LinearSystemMatrix
        if self.row_num() != self.col_num():
            return None
        n = self.row_num()
        ls = LinearSystemMatrix(self, 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._values]
        return Matrix(invA)

    def lu(self):
        # 方阵的LU分解
        assert self.row_num() == self.col_num(), "矩阵必须为方阵"
        n = self.row_num()
        A = [self.row_vector(i) for i in range(n)]
        L = [[1 if i == j else 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)])

    def rank(self):
        # 矩阵的秩
        from .LinearSystemFunc import LinearSystemFunc
        zero = Vector.zero(self.row_num())
        ls = LinearSystemFunc(self, zero)
        ls.gauss_jordan_elimination()
        zero = Vector.zero(self.col_num() + 1)
        return sum([row != zero for row in ls.Ab])

    def row_base(self):
        # 行空间的一组基
        from .LinearSystemFunc import LinearSystemFunc
        zero = Vector.zero(self.row_num())
        ls = LinearSystemFunc(self, zero)
        ls.gauss_jordan_elimination()
        zero = Vector.zero(self.col_num() + 1)
        return [new_row for new_row in [Vector(row[:self.col_num()]) for row in ls.Ab if row != zero]]

    def col_base(self):
        # 列空间的一组基
        from .LinearSystemFunc import LinearSystemFunc
        zero = Vector.zero(self.row_num())
        ls = LinearSystemFunc(self, zero)
        ls.gauss_jordan_elimination()
        return [self.col_vector(i) for i in ls.pivots]

    def dimension_null_space(self):
        # 零空间的维度
        return self.col_num() - self.rank()

    def null_base(self):
        # 零空间的一组基
        from .LinearSystemFunc import LinearSystemFunc
        zero = Vector.zero(self.row_num())
        ls = LinearSystemFunc(self, zero)
        ls.gauss_jordan_elimination()
        zero_base_list = []
        for i in range(self.col_num()):
            if i not in ls.pivots:
                temp = Vector.zero(self.col_num()).underlying_list()
                col = Matrix(ls.Ab).col_vector(i)
                for j in range(len(temp)):
                    if j < len(col) and not(is_zero(col[j])):
                        temp[j] = -col[j]
                    if i == j:
                        temp[j] = 1
                zero_base_list.append(temp)
        return [Vector(item) for item in zero_base_list]

    def dimension_left_null_space(self):
        # 左零空间的维度
        return self.row_num() - self.rank()

    def left_null_base(self):
        # 左零空间的一组基
        return self.trans().null_base()

    @classmethod
    def gram_schmidt_process(cls, basis):
        # 根据空间的一组基生成空间的一组正交基
        matrix = Matrix(basis)
        assert matrix.rank() == len(basis), \
            "列表中的向量并非空间的一组基"
        res = [basis[0]]
        for i in range(1, len(basis)):
            p = basis[i]
            for r in res:
                p = r.rectangular(p)
            res.append(p)
        return res

    @classmethod
    def arnoldi(cls, A, v, k):
        # Arnoldi算法
        u = []
        h = [[0 for col in range(k)] for row in range(k + 1)]
        w = []
        u.append(v.normalize())
        for j in range(k):
            w.insert(j, A.dot(u[j]))
            for i in range(j + 1):
                h[i][j] = u[i].dot(A.dot(u[j])) / u[i].norm()
            for i in range(j + 1):
                w[j] -= h[i][j] * u[i]
            h[j + 1][j] = w[j].norm()
            if is_zero(h[j + 1][j]):
                break
            u.append(w[j].normalize())
        return u, h

    @classmethod
    def orthonormal_basis(cls, basis):
        # 根据空间的一组基生成空间的一组标准正交基
        bas = cls.gram_schmidt_process(basis)
        return [row.normalize() for row in bas]

    @classmethod
    def orthonormal_matrix(cls, basis):
        # 由一组标准正交基生成一个标准正交矩阵
        return Matrix([[basis[i][j] for i in range(len(basis[0]))] for j in range(len(basis))])

    def qr(self):
        # 矩阵的QR分解
        assert self.row_num() == self.col_num(), \
            "矩阵必须为一个方阵"
        assert self.rank() == self.col_num(), \
            "矩阵的列向量必须线性无关"
        basis = [self.col_vector(i) for i in range(self.col_num())]
        sp = Matrix.orthonormal_basis(basis)
        Q = Matrix.orthonormal_matrix(sp)
        R = Q.trans().dot(self)
        return Q, R

    @classmethod
    def _build_col_matrix(cls, basis):
        # 根据空间的一组基建立一个列向量矩阵
        matrix = Matrix([[basis[i][j] for i in range(len(basis[0]))] for j in range(len(basis))])
        assert matrix.row_num() == matrix.col_num(), \
            "矩阵必须为一个方阵"
        assert matrix.rank() == len(basis), \
            "列表中的向量并非空间的一组基"
        return matrix

    @classmethod
    def other_to_standard(cls, other_basis, other_vector):
        # 其他坐标系坐标转换标准坐标系坐标
        matrix = cls._build_col_matrix(other_basis)
        return matrix.dot(other_vector)

    @classmethod
    def standard_to_other(cls, other_basis, vector):
        # 标准坐标系坐标转换其他坐标系坐标
        matrix = cls._build_col_matrix(other_basis)
        return matrix.inv_func().dot(vector)

    @classmethod
    def one_to_another(cls, basis, another_basis, vector):
        # 一个坐标系坐标转换另一个坐标系坐标
        standard_vector = cls.other_to_standard(basis, vector)
        return cls.standard_to_other(another_basis, standard_vector)

    def __add__(self, another):
        # 返回两个矩阵的加法结果
        assert self.shape() == another.shape(), \
            "加法错误，两个矩阵的形状必须相同"
        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(), \
            "减法错误，两个矩阵的形状必须相同"
        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 dot(self, another):
        if isinstance(another, Vector):
            # 矩阵和向量的乘法
            assert self.col_num() == len(another), \
                "矩阵列数需要与向量的长度相等"
            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(), \
                "当前矩阵列数必须等于点乘矩阵行数"
            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 __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 __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 __setitem__(self, pos, value):
        # 设置pos位置的元素
        r, c = pos
        self._values[r][c] = value

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

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

    __len__ = row_num

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

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

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

    __str__ = __repr__

from playLA.Matrix import Matrix
from playLA.Vector import Vector

if __name__ == "__main__":

    matrix = Matrix([[1, 2, 3], [3, 4, 5]])
    print(matrix)
    print("matrix.shape = {}".format(matrix.shape()))
    print("matrix.size = {}".format(matrix.size()))
    print("len(matrix) = {}".format(len(matrix)))
    print("matrix[0][0] = {}".format(matrix[0, 0]))
    print(matrix.row_vector(1))
    print(matrix.col_vector(0))
    matrix2 = Matrix([[7, 8, 9], [10, 20, 30]])
    print("add: {}".format(matrix + matrix2))
    print("subtract: {}".format(matrix2 - matrix))
    print("scalar-mul: {}".format(matrix * 3))
    print("scalar-mul: {}".format(3 * matrix))
    print("zero_2_3: {}".format(Matrix.zero(2, 3)))
    vec = Vector([23, 65, 79])
    print("matrix.dot(vec) = {}".format(matrix.dot(vec)))
    matrix3 = Matrix([[7, 9], [11, 33], [42, 97]])
    print("matrix.dot(matrix3) = {}".format(matrix.dot(matrix3)))
    print("matrix.trans = {}".format(matrix.trans()))
    I = Matrix.identity(4)
    print(I)
    idenTest = Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
    print(idenTest.dot(I))
    print(I.dot(idenTest))
    A = Matrix([[1, 2], [3, 4]])
    print(A.inv_func())
    print(A.inv_matrix())
    B = Matrix([[1, 2, 3], [4, 5, 6], [3, -3, 5]])
    L, U = B.lu()
    print(L)
    print(U)
    print(L.dot(U))
    N = Matrix([[2, 0, 0], [-1, 0, 0], [0, 0, 1]])
    print(N.rank())
    print(N.row_base())
    print(N.col_base())
    print(N.dimension_null_space())
    print(N.null_base())
    M = Matrix([[-1, 2, 3], [1, -4, -13], [-3, 5, 4]])
    print(M.rank())
    print(M.row_base())
    print(M.col_base())
    print(M.dimension_null_space())
    print(M.null_base())
    F = Matrix([[1, 2, 3, 4, 5, 6], [27, 28, 29, 30, 31, 32], [15, 16, 17, 18, 19, 20],
                [31, 32, 33, 34, 35, 36], [45, 46, 47, 48, 49, 50]])
    print(F.rank())
    print(F.row_base())
    print(F.col_base())
    print(F.dimension_null_space())
    print(F.null_base())
    print(F.dimension_left_null_space())
    print(F.left_null_base())
    basis = [Vector([2, 1]), Vector([1, 1])]
    res = Matrix.gram_schmidt_process(basis)
    print(res)
    res = Matrix.orthonormal_basis(basis)
    print(res)
    orthonormal_matrix = Matrix.orthonormal_matrix(res)
    print(orthonormal_matrix)
    print(orthonormal_matrix.dot(orthonormal_matrix.trans()))
    A1 = Matrix([[1, 1, 2], [1, 1, 0], [1, 0, 0]])
    Q, R = A1.qr()
    print(Q)
    print(R)
    print(Q.dot(R))
    other_basis = [Vector([4, 1]), Vector([2, 3])]
    other_vector = Vector([2, 2])
    vector_stard = Matrix.other_to_standard(other_basis, other_vector)
    print(vector_stard)
    print(Matrix.standard_to_other(other_basis, vector_stard))
    another_basis = [Vector([1, 3]), Vector([2, 4])]
    vector_another = Matrix.one_to_another(other_basis, another_basis, other_vector)
    print(vector_another)
    print(Matrix.other_to_standard(another_basis, vector_another))
    A2 = Matrix([[1, 2], [3, 4]])
    v = Vector([2, 2])
    u, h = Matrix.arnoldi(A2, v, 2)
    print(u)
    print(h)

运行结果

Matrix([[1, 2, 3], [3, 4, 5]])
matrix.shape = (2, 3)
matrix.size = 6
len(matrix) = 2
matrix[0][0] = 1
(3, 4, 5)
(1, 3)
add: Matrix([[8, 10, 12], [13, 24, 35]])
subtract: Matrix([[6, 6, 6], [7, 16, 25]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
scalar-mul: Matrix([[3, 6, 9], [9, 12, 15]])
zero_2_3: Matrix([[0, 0, 0], [0, 0, 0]])
matrix.dot(vec) = (390, 724)
matrix.dot(matrix3) = Matrix([[155, 366], [275, 644]])
matrix.trans = Matrix([[1, 3], [2, 4], [3, 5]])
Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
Matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])
Matrix([[-1.9999999999999996, 0.9999999999999998], [1.4999999999999998, -0.4999999999999999]])
Matrix([[-1.9999999999999996, 0.9999999999999998], [1.4999999999999998, -0.4999999999999999]])
Matrix([[1, 0, 0], [4.0, 1, 0], [3.0, 3.0, 1]])
Matrix([[1, 2, 3], [0.0, -3.0, -6.0], [0.0, 0.0, 14.0]])
Matrix([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [3.0, -3.0, 5.0]])
2
[Vector([1.0, 0.0, 0.0]), Vector([0.0, 0.0, 1.0])]
[Vector([2, -1, 0]), Vector([0, 0, 1])]
1
[Vector([0, 1, 0])]
2
[Vector([1.0, 0.0, 7.0]), Vector([-0.0, 1.0, 5.0])]
[Vector([-1, 1, -3]), Vector([2, -4, 5])]
1
[Vector([-7.0, -5.0, 1])]
2
[Vector([1.0, 0.0, -1.0000000000000007, -2.0000000000000018, -3.0000000000000013, -4.000000000000002]), Vector([0.0, 1.0, 2.0000000000000004, 3.000000000000001, 4.000000000000001, 5.000000000000002])]
[Vector([1, 27, 15, 31, 45]), Vector([2, 28, 16, 32, 46])]
4
[Vector([1.0000000000000007, -2.0000000000000004, 1, 0, 0, 0]), Vector([2.0000000000000018, -3.000000000000001, 0, 1, 0, 0]), Vector([3.0000000000000013, -4.000000000000001, 0, 0, 1, 0]), Vector([4.000000000000002, -5.000000000000002, 0, 0, 0, 1])]
3
[Vector([-0.4615384615384617, -0.5384615384615384, 1, 0, 0]), Vector([0.1538461538461533, -1.1538461538461537, 0, 1, 0]), Vector([0.6923076923076934, -1.6923076923076923, 0, 0, 1])]
[Vector([2, 1]), Vector([-0.19999999999999996, 0.4])]
[Vector([0.8944271909999159, 0.4472135954999579]), Vector([-0.44721359549995787, 0.894427190999916])]
Matrix([[0.8944271909999159, -0.44721359549995787], [0.4472135954999579, 0.894427190999916]])
Matrix([[0.9999999999999999, 0.0], [0.0, 1.0000000000000002]])
Matrix([[0.5773502691896258, 0.4082482904638628, 0.7071067811865475], [0.5773502691896258, 0.4082482904638628, -0.7071067811865476], [0.5773502691896258, -0.8164965809277263, 0.0]])
Matrix([[1.7320508075688776, 1.1547005383792517, 1.1547005383792517], [-6.661338147750939e-16, 0.8164965809277256, 0.8164965809277256], [-1.1102230246251565e-16, -1.1102230246251565e-16, 1.414213562373095]])
Matrix([[0.9999999999999999, 0.9999999999999997, 1.9999999999999996], [1.0, 0.9999999999999999, -2.220446049250313e-16], [1.0000000000000007, 3.3306690738754696e-16, 3.3306690738754696e-16]])
(12, 8)
(1.9999999999999996, 2.0)
(-15.999999999999995, 13.999999999999996)
(11.999999999999998, 8.0)
[Vector([0.7071067811865475, 0.7071067811865475]), Vector([-0.7071067811865476, 0.7071067811865474])]
[[5.0, 0.999999999999999], [1.9999999999999996, -5.551115123125784e-16], [0, 1.8054905945397061e-16]]

• 第三步，全正交化方法(FOM)

根据约束条件求解Ax=f的近似解，其中x,f为n维向量，A为n阶方阵

约束条件为：

1. 要求近似解的构成为
2. 要求在某个子空间中寻找
3. 要求找出"最好的"

注记1:x0是任给的初始解，而可以理解为"修正量"，即近似解是在初始解的基础上经过迭代修正而获得，下标k表示当前的迭代次数。修正量不是在全空间中去找，而是在的某个k维子空间中寻找，该子空间限制了寻解的范围，但它具有良好的性质与结构，因而能够期望较小的代价获得具有适当精度的近似解。当迭代次数增加时，子空间的维度也在增加，直至充满整个空间。一般希望用较小的迭代次数来获得预设的精度要求。在每一步迭代中，修正量通常使用某种最佳准则来获取。

注记2:以上的问题及约束条件规定了一大类求解方法的总框架，之下各种求解方法的区别，仅仅在于选取了何种子空间以及何种最佳准则。每一种求解方法在实现细节上的变化又派生出各种变体算法。如果把子空间取为Krylov子空间，就得到Krylov子空间方法。下面要讲的全正交化方法(FOM)以及GMRES方法都属于Krylov子空间方法，它们的区别在于使用了不同的最佳准则。

注记3:如果遇到了向量空间，首先应想到它的基，而标准正交基是性质最好的基。一旦有了基，就可以通过格拉姆-施密特方法获得标准正交基。对于Krylov子空间，当A非对称时(对称矩阵之外的矩阵，关于对称矩阵请参考线性代数整理(三) )，使用Arnoldi方法获得它们的标准正交基。不难预料，求解一般线性方程组的全正交化方法(FOM)以及GMRES方法都会用到Arnoldi方法。

设x0为任给的初始解，一般用x0替换Ax=f中的x并不能使等式成立(否则已找到解)。令，称为初始余量。由初始余量r0和系数矩阵A，可以得到k维Krylov子空间。设第k次迭代获得的解为。一般用替换Ax=f中的x并不能使等式成立(否则已找到解)，令，这是第k次迭代后的余量。现在按照约束条件1，将代入上式，有，即

            (3)

按照约束条件2，，于是它可由的基线性表示，即，于是，即，当在中发生变化时，在中发生相应的变化，从而由关系式(3)所确定的的大小和方向也发生相应的变化。很自然地，我们希望的模能够尽可能小(余量的模越小则对应的近似解越接近于精确解)，为此只需调整使得垂直于，则此时的就是中所能找到的"最好的"，即如约束条件3所要求的那样。

归结的问题：寻找，使得垂直于

由于未知，则和也是未知。这样，由条件"垂直于"仍然无法确定。这意味着该条件不够充分，注意到，则要使"垂直于"，只需令垂直于，即要求垂直于中的所有向量(包括)，于是得到如下强化的问题。

强化的问题：寻找，使得垂直于(从而使得垂直于)

现在我们来看一下如果寻找，使得垂直于会如何？因为只是中的一(大)部分，即不必垂直于中的向量。换句话说，这里要求于标准正交向量组垂直，而不要求与垂直(这样，原来的最佳准则变成了"次佳的"，见注记4)。在这样的条件下求解，就是全正交方法(FOM)所做的事情。具体求解过程如下。

既然，则可由的正交基表示，即，令，则有。现在只需确定。由垂直于，有(这里可以理解成•=0)，将(3)代入，整理得，(这里可以理解成•=•)或写成矩阵形式，将代入，得，。于是有矩阵形式。结合(2)式，得，当可逆时，就有，注意到(这里v1为r0单位化后的单位向量)，及的正交性，有，其中，。将化简的代入的表达式，得，该式意味着用乘以的第一列就得到。这样，由就得到方程组的近似解。

伪代码

ß = ||r0||
v[1] = r0 / ß      //r0单位化
for (int j = 1;j < k;j++) {    //k维子空间
    w[j] = Av[j]               //生成下一个向量
    for (int i = 1;i <= j;i++) {
        h[i][j] = (w[j],v[i])   //获取Av[j]在v[i]上的投影向量的模
    }
    for (int i = 1;i <= j;i++) {
        w[j] -= h[i][j]v[i]    //获取j+1维空间与j维体的正交向量w[j]
    }
    h[j+1][j] = ||w[j]||       //获取该正交向量w[j]的模
    if (h[j+1][j] == 0) {
        break
    }
    v[j+1] = w[j] / h[j+1][j]  //w[j]单位化
}
yk = ßh.inv()•e1     //yk=ß(Hk^-1)e1，其中e1为一维单位矩阵
xk = x0 + v•yk             //xk=x0+Vkyk

注记4:条件减弱的结果是，问题的解也被削弱，即当前条件下，不能保证与垂直，从而的一般达不到最小。实际上，由于，则在当前条件下垂直于r0(而不是).将看到，这样的是利用正交关系所能确定的中"第二最好的"或"次佳的"，而此时r0恰好是在上的投影(高维空间中的向量向低维空间投影)。换句话说，在"归结的问题"和"强化的问题"中，寻找的是"第一最好的"解，在那里是r0在中的投影(低维空间中的向量r0向高维空间投影)，而在"第二最好的"情形下，情况刚好反过来。

Java实现

public abstract class LinearSystem<T extends Number> implements FuncGroup<T> {
    //系数矩阵的行数
    protected int matrixRow;
    //系数矩阵的列数
    protected int matrixCol;
    //增广矩阵
    @Getter
    protected Ranks<T> Ab;
    //主元列
    @Getter
    protected List<Integer> pivots = new ArrayList<>();

    @SuppressWarnings("unchecked")
    public LinearSystem(Ranks<T> A, Attribute<T> b) {
        if (A.rowNum() != b.len()) {
            throw new IllegalArgumentException("系数矩阵的行数必须等于常数向量的长度");
        }
        matrixRow = A.rowNum();
        matrixCol = A.colNum();
        T[][] values = (T[][]) new Number[A.rowNum()][A.colNum() + 1];
        for (int i = 0; i < A.rowNum(); i++) {
            for (int j = 0; j <= A.colNum(); j++) {
                if (j < A.colNum()) {
                    values[i][j] = A.get(new int[]{i, j});
                } else {
                    values[i][j] = b.get(i);
                }
            }
        }
        Ab = new Matrix<>(values);
    }

    @SuppressWarnings("unchecked")
    public LinearSystem(Ranks<T> A, Ranks<T> b) {
        if (A.rowNum() != b.colNum()) {
            throw new IllegalArgumentException("系数矩阵的行数必须等于单位矩阵的列数");
        }
        matrixRow = A.rowNum();
        matrixCol = A.colNum();
        T[][] values = (T[][])new Number[A.rowNum()][A.colNum() + b.colNum()];
        for (int i = 0; i < A.rowNum(); i++) {
            for (int j = 0; j < A.colNum() + b.colNum(); j++) {
                if (j < A.colNum()) {
                    values[i][j] = A.get(new int[]{i, j});
                }else {
                    values[i][j] = b.get(new int[]{i, j - A.colNum()});
                }
            }
        }
        Ab = new Matrix<>(values);
    }

    @Override
    public boolean gaussJordanElimination(Addr<T> addr) {
        forward(addr);
        backword(addr);
        //无主元的行的非系数元素如果为0则无解
        for (int i = pivots.size(); i < matrixRow; i++) {
            if (!addr.equals(Ab.get(new int[]{i,matrixCol}),addr.zero())) {
                return false;
            }
        }
        return true;
    }

    @Override
    public void fancyPrint() {
        StringBuilder builder = new StringBuilder();
        if (Ab.colNum() - matrixCol > 1) {
            for (int i = 0; i < matrixRow; i++) {
                for (int j = 0; j < Ab.colNum(); j++) {
                    if (j < matrixCol) {
                        builder.append(" ");
                        builder.append(Ab.get(new int[]{i, j}));
                    } else if (j == matrixCol) {
                        builder.append(" | ");
                        builder.append(Ab.get(new int[]{i, j}));
                    } else {
                        builder.append(" ");
                        builder.append(Ab.get(new int[]{i, j}));
                        builder.append("\n");
                    }
                }
            }
        }else if (Ab.colNum() - matrixCol == 1) {
            for (int i = 0; i < matrixRow; i++) {
                for (int j = 0; j < Ab.colNum(); j++) {
                    if (j < matrixCol) {
                        builder.append(" ");
                        builder.append(Ab.get(new int[]{i, j}));
                    } else if (j == matrixCol) {
                        builder.append(" | ");
                        builder.append(Ab.get(new int[]{i, j}));
                        builder.append("\n");
                    }
                }
            }
        }
        System.out.println(builder.toString());
    }

    protected abstract void forward(Addr<T> addr);

    protected abstract void backword(Addr<T> addr);

    /**
     * 获取主元最大值所在的行
     * @param addr
     * @return
     */
    protected int maxRow(int indexI,int indexJ,int n,Addr<T> addr) {
        T best = Ab.get(new int[]{indexI,indexJ});
        int ret = indexI;
        for (int i = indexI + 1; i < n; i++) {
            if (addr.compair(Ab.get(new int[]{i,indexJ}),best) == 1) {
                best = Ab.get(new int[]{i,indexJ});
                ret = i;
            }
        }
        return ret;
    }

    /**
     * FOM算法
     * @param A 线性系统的系数矩阵
     * @param f 线性系统的常数向量
     * @param x0 初始解向量
     * @param k 迭代次数
     * @return 线性系统的真实解向量
     */
    public static Attribute<Double> fom(Ranks<Double> A,Attribute<Double> f,Attribute<Double> x0,int k,Addr<Double> addr) {
        //获取初始余量
        Attribute<Double> r0 = f.sub(A.dot(x0,addr),addr);
        //将初始余量放入arnoldi算法中，获取Krylov子空间的标准正交向量以及海森堡矩阵
        Tuple<List<Attribute<Double>>,Double[][]> arnoldi = Matrix.arnoldi(A,r0,k,addr);
        Double beta = r0.norm(addr);
        List<Attribute<Double>> normalBasis = arnoldi.getField(0);
        Double[][] h = arnoldi.getField(1);
        //去掉海森堡矩阵的最后一行
        Double[][] values = new Double[h.length - 1][h[0].length];
        for (int i = 0; i < values.length; i++) {
            values[i] = Arrays.copyOf(h[i],h[0].length);
        }
        Ranks<Double> H = new Matrix<>(values);
        //构建e1向量[1,0,0,0,...]
        Double[] eValue = new Double[k];
        Arrays.fill(eValue,addr.zero());
        eValue[0] = addr.one();
        Attribute<Double> e1 = new Vector<>(eValue);
        //将标准正交向量变成标准正交矩阵
        Double[][] V = new Double[normalBasis.size()][normalBasis.get(0).len()];
        for (int i = 0; i < normalBasis.size(); i++) {
            for (int j = 0; j < normalBasis.get(0).len(); j++) {
                V[i][j] = normalBasis.get(i).get(j);
            }
        }
        Ranks<Double> Vk = new Matrix<>(V);
        //获取海森堡矩阵的逆
        Ranks<Double> Hinv = H.invFunc(addr);
        if (Hinv != null) {
            //计算线性系统的真实解向量
            return x0.add(Vk.dot(Hinv.dot(e1, addr).mul(beta, addr), addr), addr);
        }
        return null;
    }

    public static void main(String[] args) {
        Ranks<Double> A = new Matrix<>(new Double[][]{{1.0,2.0},{3.0,3.0}});
        Attribute<Double> b = new Vector<>(new Double[]{3.0,7.0});
        Attribute<Double> x0 = new Vector<>(new Double[]{1.0,1.0});
        Addr<Double> addr = new Addr<Double>() {
            @Override
            public Double zero() {
                return 0.0;
            }

            @Override
            public Double one() {
                return 1.0;
            }

            @Override
            public Double add(Double lhs, Double rhs) {
                return lhs + rhs;
            }

            @Override
            public Double sub(Double lhs, Double rhs) {
                return lhs - rhs;
            }

            @Override
            public Double mul(Double k, Double hs) {
                return k * hs;
            }

            @Override
            public Double square(Double hs) {
                return hs * hs;
            }

            @Override
            public double prescription(Double hs) {
                return Math.sqrt(hs);
            }

            @Override
            public double div(Double hs, double nhs) {
                return hs / nhs;
            }

            @Override
            public Double div(Double hs, Double nhs) {
                return hs / nhs;
            }

            @Override
            public double acos(double hs) {
                return Math.acos(hs);
            }
            @Override
            public double toDegree(double hs) {
                return Math.toDegrees(hs);
            }

            @Override
            public int compair(Double lhs, Double rhs) {
                if (lhs - rhs > 0) {
                    return 1;
                }else if (lhs - rhs < 0) {
                    return -1;
                }else {
                    return 0;
                }
            }

            @Override
            public boolean equals(Double lhs, Double rhs) {
                double dis = 1e-8;
                if (Math.abs(lhs - rhs) < dis) {
                    return true;
                }
                return false;
            }

            @Override
            public Double neg(Double hs) {
                return -hs;
            }
        };
        System.out.println(LinearSystem.fom(A,b,x0,2,addr));
        FuncGroup<Double> linearSystem = new LinearSystemFunc<>(A,b);
        linearSystem.gaussJordanElimination(addr);
        linearSystem.fancyPrint();
    }
}

运行结果

Vector{values=[1.6666666666666665, 0.6666666666666667]}
 1.0 0.0 | 1.666666666666667
 0.0 1.0 | 0.6666666666666665

Python代码

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

class LinearSystemFunc:

    def __init__(self, A, b):
        assert A.row_num() == len(b), "矩阵的行数必须等于向量的长度"
        self._matrix_row = A.row_num()
        self._matrix_col = A.col_num()
        if isinstance(b, Vector):
            self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]])
                       for i in range(self._matrix_row)]
        if isinstance(b, Matrix):
            self.Ab = [Vector(A.row_vector(i).underlying_list() + b.row_vector(i).underlying_list())
                       for i in range(self._matrix_row)]
        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._matrix_row and k < self._matrix_col:
            # 看Ab[i][k]位置是否可以为主元
            max_row = self._max_row(i, k, self._matrix_row)
            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._matrix_row):
                    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][k] * self.Ab[i]

    def gauss_jordan_elimination(self):
        # 高斯-约旦消元法，如果有解，返回True;如果没有解，返回False
        self._forward()
        self._backward()
        for i in range(len(self.pivots), self._matrix_row):
            if not is_zero(self.Ab[i][-1]):
                return False
        return True

    def fancy_print(self):
        # 打印增广矩阵
        for i in range(self._matrix_row):
            print(" ".join(str(self.Ab[i][j]) for j in range(self._matrix_col)), end=" ")
            print("|", self.Ab[i][-1])

    @classmethod
    def fom(cls, A, f, x0, k):
        # FOM算法
        r0 = f - A.dot(x0)
        beta = r0.norm()
        normal_basis, h = Matrix.arnoldi(A, r0, k)
        values = [row[:] for row in h[:len(h) - 1]]
        H = Matrix(values)
        e_value = [0 for i in range(k)]
        e_value[0] = 1
        e1 = Vector(e_value)
        Vk = Matrix(normal_basis)
        H_inv = H.inv_func()
        if H_inv is not None:
            return x0 + Vk.dot(beta * H_inv.dot(e1))
        return None

from playLA.Matrix import Matrix
from playLA.Vector import Vector
from playLA.LinearSystemFunc import LinearSystemFunc

if __name__ == "__main__":

    A = Matrix([[1, 2], [3, 3]])
    b = Vector([3, 7])
    x0 = Vector([1, 1])
    print(LinearSystemFunc.fom(A, b, x0, 2))
    line_system = LinearSystemFunc(A, b)
    line_system.gauss_jordan_elimination()
    line_system.fancy_print()

运行结果

Vector{values=[1.6666666666666665, 0.6666666666666667]}
 1.0 0.0 | 1.666666666666667
 0.0 1.0 | 0.6666666666666665