Lapack和Matlab中的LU分解不同

xmakbtuz  于 2023-06-06  发布在  Matlab
关注(0)|答案(2)|浏览(177)

我正在用Lapack计算LU分解并比较结果。我使用dgetrf或dgetf2。将结果与Matlab中相应的函数lu(...)进行比较,只要输入矩阵是正方形或列数多于行数,它们就相同。对于行数多于列数的矩阵,右下方的lu矩阵的值不同。例如,Lapack中的dgetrf和dgetf2用于

A=[
    1  2  3;
    4  5  6;
    7  8  9;
   10 11 12]

给予

lu = [
   10.0 11.0    12.0; 
    0.1  0.9     1.8;
    0.7  0.33333 0.0;
    0.4  0.66666 0.45989304812834225]

但在Matlab中计算的结果相同

lu = [
   10.0 11.0    12.0; 
    0.1  0.9     1.8;
    0.7  0.33333 0.0;
    0.4  0.66666 0.0]

我不知道最后一行的最后一个值的差异来自哪里。当具有甚至更多行时,存在相同的效果。我在Swift中使用Lapack,但我不认为它与编程语言有任何关系。有线索吗?谢谢你的帮助。
我尝试改变dgetrf和dgetf2的输入参数,但效果仍然相同。

pn9klfpd

pn9klfpd1#

我写了一些MATLAB m代码来测试我的win PC运行R2021 a。我还编写了C代码来直接调用MATLAB LAPACK库的detrf()和dgetf 2()。以下是我得到的结果:

>> lu_test
A =
     1     2     3
     4     5     6
     7     8     9
    10    11    12
LU =
   10.0000   11.0000   12.0000
    0.1000    0.9000    1.8000
    0.7000    0.3333    0.0000
    0.4000    0.6667    0.4599
LUdgetrf =
   10.0000   11.0000   12.0000
    0.1000    0.9000    1.8000
    0.7000    0.3333    0.0000
    0.4000    0.6667    0.4599
LUdgetf2 =
   10.0000   11.0000   12.0000
    0.1000    0.9000    1.8000
    0.7000    0.3333    0.0000
    0.4000    0.6667    0.4599

因此,MATLAB lu()函数产生的结果与直接调用MATLAB LAPACK库函数getrf()和dgetf 2()相同,并且它们都填充了右下角的斑点。也许你可以将你的代码与下面的代码进行比较,以帮助确定你的终端出了什么问题。
m代码:

A=[ 1  2  3;
    4  5  6;
    7  8  9;
    10 11 12] 
LU = lu(A)
LUdgetrf = dgetrf(A)
LUdgetf2 = dgetf2(A)

getrf的C代码:

/* MATLAB mex routine, Calls LAPACK dgetrf routine for LU factoring     */
/* Programmer: James Tursa                                              */

/* Includes ----------------------------------------------------------- */

#include "mex.h"
#include "blas.h"
#include "lapack.h"

/* Gateway ------------------------------------------------------------ */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSignedIndex M, N, LDA, INFO, MN;
    mwSignedIndex *IPIV;
    mxArray *A;
    
    if( nlhs > 3 ) {
        mexErrMsgTxt("Too many outputs.");
    }
    if( nrhs != 1 || !mxIsDouble(prhs[0]) || mxIsSparse(prhs[0]) || mxIsComplex(prhs[0])
                  || mxGetNumberOfDimensions(prhs[0]) != 2 ) {
        mexErrMsgTxt("Need one full real double 2D matrix input.");
    }
    M = mxGetM(prhs[0]);
    N = mxGetN(prhs[0]);
    LDA = M;
    MN = M<N ? M : N;
    IPIV = mxMalloc(MN * sizeof(*IPIV));
    A = mxDuplicateArray(prhs[0]);
    dgetrf( &M, &N, mxGetData(A), &LDA, IPIV, &INFO );
    if( INFO ) {
        mexErrMsgTxt("dgetrf returned an error code.");
    }
    if( nrhs <= 1 ) {
        plhs[0] = A;
    } else if( nrhs == 2 ) {
        mexErrMsgTxt("2 outputs not yet implemented.");
    } else {
        mexErrMsgTxt("3 outputs not yet implemented.");
    }
    mxFree(IPIV);
}

dgetf 2的C代码:

/* MATLAB mex routine, Calls LAPACK dgetf2 routine for LU factoring     */
/* Programmer: James Tursa                                              */

/* Includes ----------------------------------------------------------- */

#include "mex.h"
#include "blas.h"
#include "lapack.h"

/* Gateway ------------------------------------------------------------ */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSignedIndex M, N, LDA, INFO, MN;
    mwSignedIndex *IPIV;
    mxArray *A;
    
    if( nlhs > 3 ) {
        mexErrMsgTxt("Too many outputs.");
    }
    if( nrhs != 1 || !mxIsDouble(prhs[0]) || mxIsSparse(prhs[0]) || mxIsComplex(prhs[0])
                  || mxGetNumberOfDimensions(prhs[0]) != 2 ) {
        mexErrMsgTxt("Need one full real double 2D matrix input.");
    }
    M = mxGetM(prhs[0]);
    N = mxGetN(prhs[0]);
    LDA = M;
    MN = M<N ? M : N;
    IPIV = mxMalloc(MN * sizeof(*IPIV));
    A = mxDuplicateArray(prhs[0]);
    dgetf2( &M, &N, mxGetData(A), &LDA, IPIV, &INFO );
    if( INFO ) {
        mexErrMsgTxt("dgetf2 returned an error code.");
    }
    if( nrhs <= 1 ) {
        plhs[0] = A;
    } else if( nrhs == 2 ) {
        mexErrMsgTxt("2 outputs not yet implemented.");
    } else {
        mexErrMsgTxt("3 outputs not yet implemented.");
    }
    mxFree(IPIV);
}

mex编译例程:

function compile_lapack(varargin)
libdir = 'microsoft';
lib_blas = fullfile(matlabroot,'extern','lib',computer('arch'),libdir,'libmwblas.lib');
lib_lapack = fullfile(matlabroot,'extern','lib',computer('arch'),libdir,'libmwlapack.lib');
mex(varargin{:},lib_blas,lib_lapack,'-largeArrayDims');
end

要编译mex例程:

>> compile_lapack('dgetf2.c')
Building with 'MinGW64 Compiler (C)'.
MEX completed successfully.

>> compile_lapack('dgetrf.c')
Building with 'MinGW64 Compiler (C)'.
MEX completed successfully.

更新--------------------------------------------------------
在线MATLAB R2023 a给出了以下结果:

>> LUA = lu(A)
LUA =
   10.0000   11.0000   12.0000
    0.1000    0.9000    1.8000
    0.7000    0.3333         0
    0.4000    0.6667         0
>> LUA(3,3)
ans =
   0
>> LUA(4,3)
ans =
   0

因此,在所使用的LU算法中确实出现了一些最近的变化,该变化为U(3,3)和L(4,3)点提供了精确的0。而我的PC MATLAB R2021 a给出:

>> LUA(3,3)
ans =
   6.9204e-16
>> LUA(4,3)
ans =
    0.4599

是的,这看起来很有趣,但这两个答案都是“正确的”,因为这些残差值都在简单浮点数值效应的范围内。
MATLAB Answers论坛上的反馈表明,此更改 * 可能 * 已在R2022 a版本中发生。

mwkjh3gx

mwkjh3gx2#

重要的是要注意,对于行数多于列数的矩阵,因式分解结果可能会有所不同。这是因为MATLAB和LAPACK可能使用不同的算法和策略来处理这种情况。当一个矩阵的行比列多(高矩阵),它被称为超定系统。在这种情况下,LAPACK例程可以执行 * 瘦QR* 因式分解而不是传统的LU因式分解。skinny QR 分解将矩阵分解为正交矩阵(Q)和上三角矩阵(R),使得A = QR,其中Q是MxN矩阵,R是NxN矩阵。这种因式分解通常用于求解最小二乘问题和计算超定系统的最小范数解。相比之下,MATLAB中的lu()函数通常使用部分旋转来执行标准LU因式分解,其中矩阵A被分解为下三角矩阵(L)和上三角矩阵(U),使得A = LU。
由于基础算法中的这些差异,高矩阵的因式分解结果在LAPACK和MATLAB之间可能不相同。然而,这两个因式分解仍然表示原始矩阵,并且可以用于进一步的计算或求解线性系统。如果你需要在LAPACK和MATLAB之间对行比列多的矩阵得到一致的结果,你可以考虑在这两种环境中使用相同的算法。在MATLAB中,您可以使用qr()函数和适当的输入选项执行 skinny QR 分解。

相关问题