LU分解（是/否部分主元法）+ 施密特（Schmidt）QR分解 + 吉文斯（Givens）QR分解 + Household QR分解 代码详解，可直接运行版本复测试用例（java8版）

最新推荐文章于 2023-03-26 14:52:58 发布
最新推荐文章于 2023-03-26 14:52:58 发布
阅读量5.6k 收藏 4
点赞数 2
文章标签： 矩阵 线性代数 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_38214318/article/details/112563981
版权

1）2）LU分解（是/否部分主元法）+

3）施密特（Schmidt）QR分解 +

4）吉文斯（Givens）QR分解 +

5）Household QR分解 代码详解

可直接运行版本复测试用例（java8版）

更改测试用例：更改主方法内数组数据，或直接重写I/O流

import java.math.BigDecimal;
import java.util.ArrayList;
import java.util.LinkedList;
import java.util.List;

public class Fractorization {

    public static void main(String args[]){
        run(LUFractorization(new double[][]{
                {2.0, 2.0, 2.0},
                {4.0, 7.0, 7.0},
                {6.0, 18.0, 22.0},
        }, 0));
        run(LUFractorization(new double[][]{
                {1.0, 2.0, -3.0, 4.0},
                {4.0, 8.0, 12.0, -8.0},
                {2.0, 3.0, 2.0, 1.0},
                {-3.0, -1.0, 1.0, -4.0}
        }, 1));
        run(SchmidtQR(new double[][]{
                {1, 0, -1},
                {1, 2, 1},
                {1, 1, -3},
                {0, 1, 1}
        }));
        run(GivensQR(new double[][]{
                {0, -20, -14},
                {3, 27, -4},
                {4, 11, -2}
        }));
        run(HouseholdQR(new double[][]{
                {1, 19, -34},
                {-2, -5, 20},
                {2, 8, 37}
        }));
    }

    public static void run(List<double[][]> list){
        if(list.size() == 3){
            System.out.println("The matrix L of A is");
            print(list.get(0));
            System.out.println("The matrix U of A is");
            print(list.get(1));
            System.out.println("The matrix P of A is");
            print(list.get(2));
            System.out.println();
        }
        else{
            System.out.println("The matrix Q of A is");
            print(list.get(0));
            System.out.println("The matrix R of A is");
            print(list.get(1));
            System.out.println();
        }
    }

    public static List<double[][]> LUFractorization(double[][] matrix, int model){
        List<double[][]> list = new ArrayList<>();
        int length = matrix.length;
        double[][] L = new double[length][length];
        double[][] U = new double[length][length];
        double[][] P = new double[length][length];
        int[] indexP = new int[length];

        for(int i = 0; i < length; i++)
            indexP[i] = i;
        for(int i = 0; i < length; i++){
            if(model == 1)
                exchange(matrix, L, indexP, i, findExchange(matrix, i));
            for(int j = i+1; j < length; j++)
                L[j][i] = eliminate(matrix, i, j);
        }
        for(int i = 0; i < length; i++) {
            for(int j = length-1; j >= i; j--)
                U[i][j] = matrix[i][j];
            if(model == 1)
                P[i][indexP[i]] = 1;
            else
                P[i][i] = 1;
            L[i][i] = 1;
        }
        list.add(L);
        list.add(U);
        list.add(P);
        return list;
    }

    public static List<double[][]> SchmidtQR(double[][] matrix){
        int m = matrix.length;
        int n = matrix[0].length;
        double[][] R = new double[n][n];
        List<double[][]> outList = new LinkedList<>();
        for(int i = 0; i < n; i++){
            double[][] u = new double[m][1];
            for(int j = 0; j < m; j++)
                u[j][0] = matrix[j][i];
            for(int j = 0; j < i; j++) {
                double[][] q = new double[m][1];
                for(int k = 0; k < m; k++)
                    q[k][0] = matrix[k][j];
                double r = matrixProduct(T(q), u)[0][0];
                u = matrixPlus(u, matrixNumProduct(r, q), false);
                R[j][i] = r;
            }
            double mode = 0;
            for(int j = 0; j < m; j++)
                mode += u[j][0]*u[j][0];
            double r = Math.pow(mode, 0.5);
            R[i][i] = r;
            u = matrixNumProduct(1/r, u);
            for(int j = 0; j < m; j++)
                matrix[j][i] = u[j][0];
        }
        outList.add(matrix);
        outList.add(R);
        return outList;
    }

    public static List<double[][]> GivensQR(double[][] matrix){
        int maxRow = matrix.length;
        int maxCol = matrix[0].length;
        int m = matrix.length;
        int n = matrix[0].length;
        List<double[][]> QList = new LinkedList<>();
        List<double[][]> RList = new LinkedList<>();
        List<double[][]> outList = new LinkedList<>();
        while(maxRow > 1 && maxCol > 1){
            int row = matrix.length;
            double[][] u = new double[row][1];
            for(int i = 0; i < row; i++)
                u[i][0] = matrix[i][0];
            List<double[][]> TList = new LinkedList<>();
            for(int i = 1; i < row; i++){
                if(Double.doubleToLongBits(u[i][0]) == Double.doubleToLongBits(0.00))
                    continue;
                double c = u[0][0] / Math.pow((u[0][0]*u[0][0])+(u[i][0]*u[i][0]), 0.5);
                double s = u[i][0] / Math.pow((u[0][0]*u[0][0])+(u[i][0]*u[i][0]), 0.5);
                double[][] T = eye(row);
                T[0][0] = T[i][i] = c;
                T[0][i] = s;
                T[i][0] = -s;
                u = matrixProduct(T, u);
                TList.add(T);
            }
            double[][] T = ((LinkedList<double[][]>)TList).poll();
            while(!TList.isEmpty()){
                double[][] temp = ((LinkedList<double[][]>)TList).poll();
                T = matrixProduct(temp, T);
            }
            QList.add(T);
            RList.add(matrixProduct(T, matrix));
            matrix = nextStep(matrixProduct(T, matrix));
            maxCol --;
            maxRow --;
        }
        outList.add(T(repairQList(QList, m, n)));
        outList.add(repairRList(RList, m, n));
        return outList;
    }

    public static List<double[][]> HouseholdQR(double[][] matrix){
        int maxRow = matrix.length;
        int maxCol = matrix[0].length;
        int m = matrix.length;
        int n = matrix[0].length;
        List<double[][]> QList = new LinkedList<>();
        List<double[][]> RList = new LinkedList<>();
        List<double[][]> outList = new LinkedList<>();
        while(maxRow > 1 && maxCol > 1){
            int row = matrix.length;
            double[][] u = new double[row][1];
            for(int i = 0; i < row; i++)
                u[i][0] = matrix[i][0];
            double mode = 0;
            for(int i = 0; i < row; i++)
                mode += u[i][0]*u[i][0];
            mode = Math.pow(mode, 0.5);
            u[0][0] -= mode;
            double coff = 2 / matrixProduct(T(u),u)[0][0];
            double[][] R = matrixPlus(eye(row), matrixNumProduct(coff, matrixProduct(u, T(u))), false);
            double[][] mid = matrixProduct(R, matrix);
            QList.add(R);
            RList.add(mid);
            matrix = nextStep(mid);
            maxCol --;
            maxRow --;
        }
        outList.add(T(repairQList(QList, m, n)));
        outList.add(repairRList(RList, m, n));
        return outList;
    }

    private static double[][] matrixProduct(double[][] A, double[][] B){
        int y = A.length;
        int x = B[0].length;
        double C[][] = new double[y][x];
        for (int i = 0; i < y; i++)
            for (int j = 0; j < x; j++)
                for (int k = 0; k < B.length; k++)
                    C[i][j] += A[i][k] * B[k][j];
        return C;
    }

    private static double[][] matrixPlus(double[][] A, double[][] B, boolean positive){
        int row = A.length;
        int col = A[0].length;
        double[][] C = new double[row][col];
        for(int i = 0; i < row; i++){
            for(int j = 0; j < col; j++){
                if(positive)
                    C[i][j] = A[i][j] + B[i][j];
                else
                    C[i][j] = A[i][j] - B[i][j];
            }
        }
        return C;
    }

    private static double[][] matrixNumProduct(double b, double[][] A){
        int row = A.length;
        int col = A[0].length;
        double[][] C = new double[row][col];
        for(int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                    C[i][j] = b * A[i][j];
            }
        }
        return C;
    }

    private static double[][] repairQList(List<double[][]> QList, int m, int n){
        double[][] Q = ((LinkedList<double[][]>)QList).poll();
        while(!QList.isEmpty()){
            double[][] temp = ((LinkedList<double[][]>)QList).poll();
            int row = temp.length;
            int col = temp[0].length;
            double[][] productor = new double[m][n];
            for(int i = 0; i < row; i++)
                for (int j = 0; j < col; j++)
                    productor[m-row+i][n-col+j] = temp[i][j];
            for(int i = 0; i < m-row; i++)
                productor[i][i] = 1;
            Q = matrixProduct(productor, Q);
        }
        return Q;
    }

    private static double[][] repairRList(List<double[][]> RList, int m, int n){
        double[][] R = ((LinkedList<double[][]>)RList).poll();
        while(!RList.isEmpty()){
            double[][] temp = ((LinkedList<double[][]>)RList).poll();
            int row = temp.length;
            int col = temp[0].length;
            for(int i = 0; i < row; i++)
                for (int j = 0; j < col; j++)
                    R[m-row+i][n-col+j] = temp[i][j];
        }
        return R;
    }

    private static double[][] T(double[][] A){
        int row = A.length;
        int col = A[0].length;
        double[][] C = new double[col][row];
        for(int i = 0; i < col; i++) {
            for (int j = 0; j < row; j++) {
                C[i][j] = A[j][i];
            }
        }
        return C;
    }

    private static double[][] eye(int x){
        double[][] I = new double[x][x];
        for(int i = 0; i < x; i++)
            I[i][i] = 1;
        return I;
    }

    private static double[][] nextStep(double[][] A){
        int row = A.length;
        int col = A[0].length;
        double[][] C = new double[row-1][col-1];
        for(int i = 0; i < col-1; i++)
            for (int j = 0; j < row-1; j++)
                C[i][j] = A[i+1][j+1];
        return C;
    }

    private static int findExchange(double[][] matrix, int pivot){
        int length = matrix.length;
        double max = matrix[pivot][pivot];
        int out = pivot;
        for(int i = pivot; i < length; i++){
            if(absoluteValue(matrix[i][pivot]) > absoluteValue(max)){
                max = matrix[i][pivot];
                out = i;
            }
        }
        return out;
    }

    private static void exchange(double[][] matrix, double[][] L, int[]indexP, int i, int exchange){
        double[] temp = matrix[i];
        matrix[i] = matrix[exchange];
        matrix[exchange] = temp;
        double[] LTemp = L[i];
        L[i] = L[exchange];
        L[exchange] = LTemp;
        int tempIndex = indexP[i];
        indexP[i] = indexP[exchange];
        indexP[exchange] = tempIndex;
    }

    private static double eliminate(double[][] matrix, int i, int j){
        int length = matrix.length;
        double coefficient = matrix[j][i]/matrix[i][i];
        matrix[j][i] = 0;
        for(int k = i+1; k < length; k++)
            matrix[j][k] -= coefficient*matrix[i][k];
        return coefficient;
    }

    private static double absoluteValue(double x){
        return x >= 0 ? x : -x;
    }

    private static double[][] copy(double[][] matrix){
        int length = matrix.length;
        double[][] out = new double[length][length];
        for(int i = 0; i < length; i++)
            for(int j = 0; j < length; j++)
                out[i][j] = matrix[i][j];
        return out;
    }


    private static void print(double[][] matrix){
        int row = matrix.length;
        int col = matrix[0].length;
        for(int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                BigDecimal bigDecimal = new BigDecimal(matrix[i][j]);
                String out = String.valueOf(bigDecimal.setScale(3, BigDecimal.ROUND_HALF_UP).doubleValue());
                int add = 8 - out.length();
                for(int k = 0; k < add; k++)
                    out += " ";
                System.out.print(out);
            }
            System.out.println();
        }
    }
}

 

pw_hu
关注
  • 2
    点赞
  • 踩
  • 4
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
Givens变换的matlab程序
05-15
数值分析中的一个重要变化GIVENS 变换，该资源包内包含详细的程序和注释
Givens变换
ddzai7557的博客
09-17 1万+
摘抄自《数值线性代数(徐树方)》 Givens变换是一种将n维向量x在第(i，k)两个维度确定的坐标平面内进行旋转（从而将其中一个分量化0）的变换，因此它又叫平面旋转变换。 与Householder变换将一个向量中的若干个分量化0相比，Givens变换将向量的其中一个分量化0。 本篇先介绍Givens变换的定义及其性质，再推导一种用于求Givens变换的数值化方法 一、Given...
吉文斯旋转
weixin_30662539的博客
11-13 798
在数值线性代数中，吉文斯旋转（Givens rotation）是在两个坐标轴所展开的平面中的旋转。吉文斯旋转得名于华莱士·吉文斯，他在 1950 年代工作于阿贡国家实验室时把它介入到数值分析中。 目录 1矩阵表示2稳定计算3参见4引用 矩阵表示 吉文斯旋转表示为如下形式的矩阵 这里的 c = cos(θ) 和 s = sin(θ) 出现在第 i 行和第...
矩阵分析与应用(17)
qq_40206924的博客
07-29 2557
学习来源《矩阵分析与应用》张贤达清华大学出版社。
线性代数 --- Gauss消元的部分主元法和完全主元法
松下J27录放机
05-18 6652
Gauss消元的部分主元法和完全主元法
部分主元法的LU分解
12-26
用MATLAB实现部分主元法的LU分解，选取列中绝对值最大的行实现行交换。
LU.zip_article22t_foundatv_wheel8em_主元的lu分解_列主元LU分解
09-23
实现LU分解，选列主元的LU分解和不选主元的LU分解，在残差允许的范围内与MATLAB内置可以达到相同效果
teacher_fortran.rar_LU_LU fortran_lu 矩阵_lu分解法_teacher
09-21
LU分解法 修正的LU分解法 高斯全选主元法 雅克比矩阵等算法
LU分解法、列主元高斯法、Jacobi迭代法、Gauss-Seidel法的原理及Matlab程序
05-09
LU分解法、列主元高斯法、Jacobi迭代法、Gauss-Seidel法的原理及Matlab程序
LU分解法matlab脚本文件
05-04
我上传的是一个matlab脚本文件，程序会自动生成一个元素大小在3到5之间的四阶非奇异矩阵，向量b也是自动生成的，元素大小在4到6之间，有需要的朋友可根据实际情况进行更改
Givens代码及说明.rar_Givens变换_givens变换matlab_givens变换代码_matlabgivens函
07-14
数值线性代数中的givens变换，matlab语言实现
三种消元法（全主元、Gauss消去法、列主元）
04-23
三种消元法（全主元、Gauss消去法、列主元）三种消元法（全主元、Gauss消去法、列主元）三种消元法（全主元、Gauss消去法、列主元）三种消元法（全主元、Gauss消去法、列主元）
矩阵分析与计算学习记录-矩阵分解
Ricardo_PING的Blog
10-27 3187
满秩分解：存在性、方法 三角分解：Doolittle分解、两种求解方法、cholesky分解 QR分解：定义、Householder变换、Givens变换、Schmidt正交化方法求QR分解、上Hessenberg矩阵 奇异值分解
线性代数 --- Gauss消元的部分主元法和完全主元法(补充)
松下J27录放机
12-22 1838
直接写出部分主元高斯消元法中的L矩阵和P矩阵的方法。
[矩阵的QR分解系列二] 吉文斯(Givens)变换
热门推荐
honyniu的专栏
11-30 2万+
吉文斯变换简介 之前介绍的矩阵的三角分解系列介绍了利用矩阵初等变换解决了矩阵三角化问题以及具体的三角分解。但是以初等变换工具的三角分解方法并不能消除病态线性方程组不稳定问题，而且有时候对于可逆矩阵有可能也不存在三角分解。所以后面为了解决这里问题，发展出来了以正交(酉)变换的矩阵的QR(正交三角)分解，矩阵的正交三角分解是一种对任何可逆矩阵均存在理想分解。进行QR分解需要用到施密特(Schmidt)正交规范化，吉文斯(Givens)变换和豪斯霍尔德(Householder)变换等。这里矩阵的QR分解系列教程主
基于Givens矩阵的QR矩阵分解
ONLY_ME_CANTHIS的博客
03-26 2812
Givens矩阵是一种旋转矩阵，可以将一个向量旋转到另一个向量的方向。在QR分解中，我们使用Givens矩阵将矩阵的列向量逐个旋转，使得矩阵变为上三角矩阵。QR分解是一种将矩阵分解为正交矩阵和上三角矩阵的方法。在QR分解中，正交矩阵Q的转置是它的逆矩阵，因此QR分解可以用于求解线性方程组、最小二乘问题等。对矩阵A1的第二列进行Givens变换，使得A1的第二列的下面的元素都变为0。对矩阵A的第一列进行Givens变换，使得A的第一列的下面的元素都变为0。这样得到一个新的矩阵A1和一个Givens矩阵G1。
QR分解－givens旋转与Householder变换
m0_37604894的博客
03-08 1万+
目录一、Householder QR分解Householder矩阵Householder变换Householder QR分解Householder物理意义二、Givens 旋转QR分解三、QR分解适用情况 一、Householder QR分解 Householder矩阵 Householder变换 同理 例子：将列向量其他元素置零 Householder QR分解 按列消元，一列一列操作 详细消元顺序见：householder进行QR分解 matlab代码： Householder物理意义 矩
部分选主元的LU分解，处理较为巧妙，性能高效
摄心神，致虚极，守静笃
09-01 4644
http://www.johnloomis.org/ece538/notes/Matrix/ludcmp.html LU Decomposition Download: matrix.zip. See: ludcmp.cpp  lusolve.cpp (test program). LU decomposition is a method of solving linear
LU分解（LU Factorization）计算方法（手算+MATLAB），关于置换矩阵（Permutation Matrix），部分主元消去法（Partial Pivoting）
Zijie123pea的博客
02-15 1万+
背景：   求解一些列具有相同系数矩阵的线性方程，如：Ax=b1Ax=b_1Ax=b1​，Ax=b2Ax=b_2Ax=b2​，…Ax=bpAx=b_pAx=bp​等。当矩阵AAA可逆时，可以先求出矩阵AAA的逆A−1A^{-1}A−1，再计算A−1b1A^{-1}b_1A−1b1​，A−1b2A^{-1}b_2A−1b2​等。   但是，也可以用LU分解法来解这一系列方程：先使用初等行变换化简解出Ax=b1Ax=b_1Ax=b1​，并同时得到矩阵AAA的LU分解，剩下的方程使用LU分解法求解即可。 方法：
用代码实现主元Lu分解
最新发布
05-16
主元Lu分解是一种矩阵分解方法，可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积，其中下三角矩阵和上三角矩阵都是对角线元素为1的三角矩阵。下面是用Python实现主元Lu分解的代码： ```python import numpy as np def lu_factor(A): n = A.shape[0] L = np.eye(n) U = np.copy(A) for k in range(n-1): pivot = np.argmax(np.abs(U[k:, k])) + k if U[pivot, k] == 0: raise ValueError("LU factorization failed: singular matrix") if pivot != k: U[[pivot, k], k:] = U[[k, pivot], k:] L[[pivot, k], :k] = L[[k, pivot], :k] L[k+1:, k] = U[k+1:, k] / U[k, k] U[k+1:, k:] -= L[k+1:, k][:, np.newaxis] * U[k, k:] return L, U ``` 这个函数接受一个$n\times n$的矩阵$A$作为输入，返回一个形如$(L,U)$的元组，其中$L$是下三角矩阵，$U$是上三角矩阵。这里我们采用了部分选主元的方法，即在第$k$列中选取绝对值最大的元素作为主元。如果主元为0，则无法进行分解，抛出异常。选定主元后，如果不是第$k$行，则将第$k$行和主元所在的行互换，以确保主元在第$k$行。接着，将$L$的第$k+1$行到第$n$行第$k$列的元素设置为$U$相应元素除以主元，即$L_{i,k}=U_{i,k}/U_{k,k}$。然后，将$U$的第$k+1$行到第$n$行第$k$列的元素减去$L$对应元素乘以$U$的第$k$行第$k$列的元素，即$U_{i,j} \leftarrow U_{i,j}-L_{i,k}U_{k,j}$。最后，返回$L$和$U$。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值