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 ?

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的头文件:

// 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》论文笔记

sklearn.decomposition.PCA 的简单特征向量图

交叉分解+Cross decomposition

[简单路径] Useful Decomposition

LU Decomposition