PLU Decomposition
Posted yfceshi
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了PLU Decomposition相关的知识,希望对你有一定的参考价值。
PLU分解的优点是,能够将Ax=b的矩阵,转换成Ly=b, Ux = y
的形式。当我们改变系数矩阵b时,此时因为矩阵L和U均是固定
的,所以总能高效的求出矩阵的解。
// LU.cpp : Defines the entry point for the console application. // /************************************************ * Author: JohnsonDu * From: Institute of Computing Technology * University of Chinese Academy of Science * Time: 2014-10-7 * Content: PLU decomposition *************************************************/ #include "stdafx.h" #define MAXN 5005 #define eps 1e-9 // 精度 int n, m; double mat[MAXN][MAXN]; // 输入矩阵 double matL[MAXN][MAXN]; // 矩阵L double matU[MAXN][MAXN]; // 矩阵U int matP[MAXN][MAXN]; // 矩阵P int seq[MAXN]; // 记录行变换 //double vecB[MAXN]; //double vecY[MAXN]; //double vecX[MAXN]; //double matPb[MAXN]; void menu() { printf("----------------PLU Factorization---------------\n"); printf("| Please follow the instruction |\n"); printf("| to determine the LU decomposition |\n"); printf("| PA = LU |\n"); printf("------------------------------------------------\n\n"); } void initLMatrix() { memset(matU, 0, sizeof(matU)); memset(matL, 0, sizeof(matL)); memset(matP, 0, sizeof(matP)); } void padLMatrix() { for(int i = 0; i < n; i ++) matL[i][i] = 1.0; } inline double Abs(double x) { return x < 0 ? -x : x; } void displayLU() { // 输出矩阵L printf("\n----------------------\n"); printf("Matrix L follows: \n"); for(int i = 0; i < n; i ++) { for(int j = 0; j < (n < m ?当中stdafx.h的头文件:n : m); j ++) printf("%.3f ", matL[i][j]); printf("\n"); } // 输出矩阵U printf("\nMatrix U follows: \n"); for(int i = 0; i < (n < m ? n : m); i ++) { for(int j = 0; j < m; j ++) printf("%.3f ", matU[i][j]); printf("\n"); } // 输出矩阵P printf("\nMatrix P follows: \n"); for(int i = 0; i < n; i ++) { for(int j = 0; j < n; j ++) printf("%d ", matP[i][j]); printf("\n"); } printf("----------------------\n"); } /* // 输出LU的过程及终于解 void displaySolution() { // 输出矩阵Pb printf("\nMatrix Pb follows: \n"); for(int i = 0; i < n; i ++) { printf("%.3f\n", matPb[i]); } // 输出向量y printf("\nVector Y follows: \n"); for(int i = 0; i < n; i ++) { printf("%.3f\n", vecY[i]); } printf("\n"); // 输出解向量x printf("\Vector X follows: \n"); for(int i = 0; i < n; i ++) { printf("%.3f\n", vecX[i]); } printf("\n"); } */ // 交换元素 inline void swap(int &a, int &b) { int t = a; a = b; b = t; } // 高斯消元部分 void gauss() { int i; int col; int max_r; col = 0; //处理的当前列 // 从第一行開始进行消元 // k为处理的当前行 for(int k = 0; k < n && col < min(n, m); k ++, col ++) { // 寻找当前col列的绝对值最大值 max_r = k; for(i = k + 1; i < n; i ++) if(Abs(mat[i][col]) > Abs(mat[max_r][col])) max_r = i; // 进行行交换 if(max_r != k) { for(int j = col; j < m; j ++) swap(mat[k][j], mat[max_r][j]); swap(seq[k], seq[max_r]); for(int j = 0; j < n; j ++) swap(matL[k][j], matL[max_r][j]); } // 当前主元为零, 继续 if(Abs(mat[k][col]) < eps){ continue; } // 消元部分,并获得L矩阵 for(int i = k + 1; i < n; i ++) { double t = mat[i][col] / mat[k][col]; matL[i][col] = t; for(int j = col; j < m; j ++) mat[i][j] -= t * mat[k][j]; } } // 为矩阵U进行赋值 for(int i = 0; i < n; i ++) for(int j = 0; j < m; j ++) matU[i][j] = mat[i][j]; // 生成矩阵P for(int i = 0; i < n; i ++) matP[i][seq[i]] = 1.0; // 为矩阵L加入对角线元素 padLMatrix(); } /* // 计算Pb的值 void calcPb() { for(int i = 0; i < n; i ++) matPb[i] = 0.0; //cout << "-----------" << endl; for(int i = 0; i < n; i ++) { double t = 0.0; for(int j = 0; j < n; j ++) { t = t + 1.0 * matP[i][j] * vecB[j]; //cout << t << endl; //cout << matP[i][j] * vecB[j] << "---" << endl; } matPb[i] = t; //cout << matPb[i] << endl; } } // 计算Ly = Pb, y向量 void calcY() { vecY[0] = matPb[0]; for(int i = 1; i < n; i ++) { double t = 0.0; for(int j = 0; j < i; j ++) t += vecY[j] * matL[i][j]; vecY[i] = matPb[i] - t; } } // 计算Ux = y, y向量 void calcX() { vecX[n-1] = vecY[n-1] / matU[n-1][n-1]; for(int i = n-2; i >= 0; i --) { double t = 0.0; for(int j = n-1; j > i; j --) t += vecX[j] * matU[i][j]; vecX[i] = (vecY[i] - t) / matU[i][i]; } } */ int _tmain(int argc, _TCHAR* argv[]) { menu(); while(true) { printf("Please input the matrix's dimension n & m: "); // 输入矩阵的行n和列m scanf("%d%d", &n, &m); printf("Please input the matrix: \n"); // 输入矩阵 for(int i = 0; i < n; i ++) { for(int j = 0; j < m; j ++) cin >> mat[i][j]; seq[i] = i; } // 初始化为0 initLMatrix(); // 高斯消元 gauss(); // 输出P, L, U矩阵 displayLU(); system("pause"); system("cls"); menu(); /* //此处是输入b,求取x, y 和 pb while(true){ printf("please input vector b(whose length equals to %d): \n", n); for(int i = 0; i < n; i ++) cin >> vecB[i]; calcPb(); calcY(); calcX(); displaySolution(); } */ } return 0; }
// stdafx.h : include file for standard system include files, // or project specific include files that are used frequently, but // are changed infrequently // #pragma once #define _CRT_SECURE_NO_WARNINGS #include "targetver.h" #include <stdio.h> #include <tchar.h> #include <iostream> using namespace std;
以上是关于PLU Decomposition的主要内容,如果未能解决你的问题,请参考以下文章
sklearn.decomposition.DicitonaryLearning.fit 中的 y 参数有啥作用?
《Single Image Depth Prediction with Wavelet Decomposition》论文笔记