精华内容
下载资源
问答
  • LDLT分解

    千次阅读 2018-06-07 09:59:39
    LDLTLDLTLDL^T分解 首先是一个矩阵是可以被下三角分解,即LDU A=BTMBA=BTMBA = B^T M B 判断一个矩阵是不是正定的,可以通过主行列式的值进行判断。 正定的判断是可以通过,M 的 trace>0,那么A就是正定...

    LDLT L D L T 分解
    首先是一个矩阵是可以被下三角分解,即LDU
    A=BTMB A = B T M B
    判断一个矩阵是不是正定的,可以通过主行列式的值进行判断。
    正定的判断是可以通过,M 的 trace>0,那么A就是正定的,因为可以通过M判断A是不是正定的。
    [xyz]123xyz=x2+2y2+3z2>=0 [ x y z ] [ 1 2 3 ] [ x y z ] = x 2 + 2 y 2 + 3 z 2 >= 0 这个就是判断正定可以使用的方法
    由于此时B不会是 [xyz] [ x y z ] 的形式;因此在这里,我们构造一个 [xyz]113/2111 [ x y z ] [ 1 1 1 3 / 2 1 1 ] 得到 [xx+y3x/2+y+z] [ x x + y 3 x / 2 + y + z ] 的形式,所以得到 XTAX=XTBTMBX=yTMy>=0 X T A X = X T B T M B X = y T M y >= 0 的形式,因此可以得到分解之后的操作仍然是正定的。
    所以 LDLT L D L T 分解的形式就如上所示,可以用于metric learning中间的 transform 的过程。

    展开全文
  • LDLT分解协处理器的并行结构研究.pdf
  • LU分解,LDLT分解,Cholesky分解

    千次阅读 2018-10-27 01:07:19
    LU分解 如果方阵是非奇异的,即的行列式不为0,LU分解总是存在的。 A=LU,将系数矩阵A转变成等价的两个矩阵L和U的乘积,其中L和U分别是下三角和上三角矩阵,而且要求L的对角元素都是1,形式如下: ...LDLT分解(LU...

    LU分解

    如果方阵是非奇异的,即的行列式不为0,LU分解总是存在的

    A=LU,将系数矩阵A转变成等价的两个矩阵L和U的乘积,其中L和U分别是下三角和上三角矩阵,而且要求L的对角元素都是1,形式如下:

    本质上,LU分解是高斯消元的一种表达方式。首先,对矩阵A通过初等行变换将其变为一个上三角矩阵,然后,将原始矩阵A变为上三角矩阵的过程,对应的变换矩阵为一个下三角矩阵。

    LDLT分解(LU的进一步分解)

    A为对称矩阵,那么会产生A=LDLT分解

    定理:若对称矩阵A的各阶顺序主子式不为零时,则A可以唯一分解为A= LDLT

    证:当矩阵A的各阶顺序主子式不为零时,A有唯一的Doolittle分解A= LU,矩阵U的对角线元素Uii 不等于0,将矩阵U的每行依次提出

    A= LDLT

    Cholesky分解

    如果A是正定矩阵,那么A可以唯一分解为A=LL^T,

    证:如果A是正定矩阵,那么A是对称的,且顺序主子式大于0,则可以唯一分解为A= LDLT

    将D分解为\sqrt{D}\sqrt{D}

    A=LDL^T=L\sqrt{D}\sqrt{D}L^T=L\sqrt{D}\sqrt{D}^TL^T=L\sqrt{D}(L\sqrt{D})^T,且分解唯一。

    如果A是半正定的,也可以分解,不过这时候L就不唯一了.

    参考:https://blog.csdn.net/zhouliyang1990/article/details/21952485

    展开全文
  • 解线性方程组的LDLT分解法 通过本实验加深对LDLT分解法的构造过程的理解; 2.能对上述LDLT分解法提出正确的算法描述编程实现。
  • 最近碰到求解线性方程组以及求矩阵的特征...这篇博文将给出一个对称矩阵的RtDR分解方法(书里面一般都是LDLt分解,我直接求得是转置,即R=L')。 数值计算不同于基本数学理论,《线性代数》以及《高等数学基础》是...

    最近碰到求解线性方程组以及求矩阵的特征值等问题,OpenCV自带的算法实在是太慢了,另外我还试了Eigen库,比OpenCv虽然快了一倍,但是比Matlab还是慢了一个量级不止。。。因此我决定自己编写几个程序以满足我的特定需要。这篇博文将给出一个对称矩阵的RtDR分解方法(书里面一般都是LDLt分解,我直接求得是转置,即R=L')。

    数值计算不同于基本数学理论,《线性代数》以及《高等数学基础》是工科的数学课程,里面介绍了很多一般线性代数问题的解法,但是那些只是在理论上可行,并没有考虑计算机存储的舍入误差的影响,如果照着课本上的思路来实现算法,则在计算效率和计算结果的精度上都与Matlab相差甚远。为此,我花了一些时间学习了两本数值分析的书,其中《Matrix Computations》(4th Edition)是非常经典的书了,写得也很赞。根据书中介绍的方法写了一些函数,在处理中小型矩阵(e.g.N在200以内)时,具有较快的速度和较高的精度。

    本文将给出这样一个函数:输入一个实对称矩阵A,计算它的RtDR分解:A=R'*D*R,其中R是单位上三角矩阵,D是对角矩阵。该方法要求A的所有顺序主子式都是非奇异的。

    这是Matlab代码 version#1:

    function  [L,D] = LDLtDecomp(A)
    % A = L*D*L'
    % A is symmetric
    % L is unit lower triangular
    % D is diagonal
    n = size(A,1);
    v = zeros(n,1);
    A(2:n,1) = A(2:n,1) / A(1,1);
    for i = 2:n
        v(1:i-1) = A(i,1:i-1) .* diag(A(1:i-1,1:i-1))';
        A(i,i) = A(i,i) - A(i,1:i-1) * v(1:i-1);
        A(i+1:n,i) = (A(i+1:n,i) - A(i+1:n,1:i-1) * v(1:i-1)) / A(i,i);
    end
    D = diag(diag(A));
    L = A;
    for i = 1:n
        L(i,i) = 1;
        for j = i+1:n
            L(i,j) = 0;
        end
    end

    这是Matlab代码version#2:

    function  [R,D] = RtDRdecomp(A)
    % The transposed version of LDLtDecomp()
    % A = R'*D*R
    % A is symmetric
    % R is unit upper triangular
    % D is diagonal
    assert(norm(A-A',1) < 1e-10); % A should be symmetric
    n = size(A,1);
    A(1,2:n) = A(1,2:n) / A(1,1);
    for i = 2:n
        v = A(1:i-1,i) .* diag(A(1:i-1,1:i-1));
        A(i,i:n) = A(i,i:n) - v' * A(1:i-1,i:n);
        A(i,i+1:n) = A(i,i+1:n) / A(i,i);
    end
    D = diag(diag(A));
    R = A;
    for i = 1:n
        R(i,i) = 1;
        R(i,1:i-1) = 0;
    end

    这是C代码:

    //Here we introduce another symmetric matrix decomposition that does not require a definiteness property.
    //by: yuxianguo, 2018/11/30
    int //return 0 means failure; otherwise success (require A to have an LU factorization)
    yuRtDRdecomp( //Compute the decomposion: A=R'*D*R, where R is unit upper triangular and D is diagonal; D is stored in A.diagonal, R is stored in A.upperTriangular
    	double *A, int N) //A[N*(N+1)/2] is the upper triangular part of a symmetric matrix; require all the principle submatrices to be non-singular, or the decomposition fails
    {
    	int i, j, k, Ni = N;
    	double a, b, *p = A, *q;
    	for(i = 0; i < N; i++, Ni--) {
    		//Ni = N - i; p points to A[i][i];
    		//for(j=0;j<i;j++)//here we use i-j as the loop index
    		for(j = i, q = A; j; j--) {
    			a = *q; q += j; b = *q++;
    			a *= -b; *p += a * b;
    			for(k = 1; k < Ni; k++)
    				p[k] += a * *q++;
    		}
    		a = *p++;
    		if(!a)
    			return 0; //fail
    		a = 1.0 / a;
    		for(k = 1; k < Ni; k++)
    			*p++ *= a;
    	}
    	return 1;
    }

    关于C代码的补充说明:

    (1)在我的设置里,所有对称矩阵只保存其上三角部分,一个N-by-N大小的矩阵,只需保存N*(N+1)/2个元素;

    (2)上述代码将输入矩阵A的分解结果D和R分别覆写到A的内存中,其中D为A的对角线部分,R为A的上三角部分(R对角线元素全是1);

    (3)我以前写程序喜欢用模板、SSE(AVX)、OpenMP(多核并行)、pthread(多线程)。现在更喜欢写纯C代码,因为好移植,也好修改;

    (4)上述算法不太适合处理大型矩阵,比如N>300的情形。书上介绍说处理大型矩阵一般都考虑并行算法,并行算法不仅仅是使用多核或多线程,还需要从算法设计上考虑并行,在矩阵计算方面,一般使用矩阵分块方法。我目前还没有处理大矩阵的需求,即使是大数据,往往也只有小矩阵(e.g.协方差矩阵);

    (5)我喜欢用返回值0表示函数失败,这与很多经典的逻辑不一致,因为我有很多函数都是返回的指针,在那些函数里我可以用返回NULL表示函数失败。

    下面是一个简单的例子:

    (1)用Matlab产生一个对称矩阵

    n = randi(90) + 15;
    U = rand(n) * - 0.5;
    D = rand(n,1) * 2 - 0.5;
    A = U * diag(D) * U';
    yusave('A',A);

    其中yusave是我写的一个函数,将Matlab矩阵保存到二进制文件中。

    (2)使用C代码处理上面生成的矩阵,并保存结果

    int main() {
    	Buffer Abuf; reserve(&Abuf, 0);
    	int N;
    	yuLoad("A", 1, &Abuf, &N, 0, 0, 0);
    	double *A1 = (double*)Abuf.p;
    	double *A = new double[N*(N + 1) / 2], *p = A;
    	for(int i = 0; i < N; i++)
    		for(int j = i; j < N; j++)
    			*p++ = A1[i * N + j];
    	
    	if(!yuRtDRdecomp(A, N)) {
    		printf("failed!");
    		return getchar();
    	}
    	p = A;
    	memset(A1, 0, N * N << DBLShift);
    	for(int i = 0; i < N; i++)
    		for(int j = i; j < N; j++)
    			A1[i * N + j] = *p++;
    	yuSave("X", A1, N, N, 1, DOUBLE64);
    
    	release(&Abuf);
    	delete[] A;
    	return 0;
    }

    main函数先把数据文件读进来,然后提取矩阵的上三角部分,再调用函数进行矩阵分解,最后将分解后的矩阵恢复到方阵形式以便保存到文本中。

    (3)在Matlab中查看结果

    X = ld('X');
    D = diag(diag(X));
    R = X;
    for i = 1:n
        R(i,i) = 1;
        R(i,1:i-1) = 0;
    end
    E = abs(A - R'*D*R);
    max(E(:))

    在我们的测试中,n=104,max(E(:))=5.9508e-14

    注:double数据的舍入误差大概是DBL_EPSILON=2e-16,基于double类型的编程求解算法结果精度不会比2e-16更高。另外Matlab自己也有计算误差,在上面的例子中,即使我们的R和D是完全正确的,max(E(:))也可能不是0.

     

    总结:

    最后小结一下。这篇博文旨在抛砖引玉,大家在编写数值计算程序时,最好能了解一点相关的知识。数值计算方法和传统的数学方法并不完全一样。(PS:本文的例子可能不够好,或许用QR分解作为案例来讲解会更恰当一些~~~。管他呢,代码请随便拿去用,只求别抹掉我的署名!)

    另外啰嗦一下,个人以为RtDR分解(或LDLt分解)要比Cholesky分解(A=L'*L,其中L为下三角矩阵)更有用:首先二者的计算量是相同的,都是O(n^3/3)的水平;RtDR分解只要求A对称且存在LU分解,而Cholesky分解要求A为对称正定矩阵,因此存在Cholesky分解的矩阵必然能进行RtDR分解,反之则不然;最后,RtDR分解能得到一个对角阵,由此很容易得到A的正、负特征值的个数(知识回顾:两个二次型之间存在可逆线性变换的充要条件是它们的正、负惯性指数相同)。

    展开全文
  • 2. 用LDLT分解求解方程组: x1+2x2+3x3 = -3 2x1+x2-2x3 = 10 3x1-2x2+x3 = 7
  • 第三课 LDLt分解高斯消元 仔细观察上一回的内容以及代码,我们会发现在对系数矩阵分解时,需要占用很长的时间。 如果系数矩阵为一个对称矩阵,也就是 比如 按照上一节中,可以把[A]分解为 如果把[U]中每一行都...

    第三课 LDLt分解高斯消元

    仔细观察上一回的内容以及代码,我们会发现在对系数矩阵分解时,需要占用很长的时间。
    如果系数矩阵为一个对称矩阵,也就是
    在这里插入图片描述
    比如
    在这里插入图片描述
    按照上一节中,可以把[A]分解为
    在这里插入图片描述
    在这里插入图片描述
    如果把[U]中每一行都除以主对角线的数,那么将变成在这里插入图片描述
    其实可以很容易看出
    在这里插入图片描述
    所以上三角矩阵可以表示成下列形式
    在这里插入图片描述
    其中[D]为系数矩阵主对角线构成的矩阵
    在这里插入图片描述
    因此,如果[A]是一个对称矩阵,我们能写成
    在这里插入图片描述
    [Lt]是[L]的转置
    计算实例:
    在这里插入图片描述
    [L][D][Lt]能被写成下面的形式
    在这里插入图片描述
    因此
    在这里插入图片描述
    各个系数的求解结果为在这里插入图片描述
    因此
    在这里插入图片描述
    在这里插入图片描述
    然后用‘向前迭代法’
    在这里插入图片描述
    得到
    在这里插入图片描述
    然后用‘向后迭代法’求解
    在这里插入图片描述
    得到
    在这里插入图片描述
    代码如下:还是包括一个主程序和三个子程序,子程序分别为ldlt分解的ldlt、向前迭代的ldlfor、向后迭代的subbac
    主程序

    #使用LDLT分解的高斯消元法
    import numpy as np
    import math
    import B
    n=3
    d=np.zeros((n,1))
    a=np.array([[3,-2,1],[-2,3,2],[1,2,2]],dtype=np.float)
    b=np.array([[3],[-3],[2]],dtype=np.float)
    print('系数矩阵')
    for i in range(1,n+1):
        print(a[i-1,:])
    print('右手边向量',b)
    B.ldlt(a,d)
    print('下三角')
    for i in range(1,n+1):
        print(a[i-1,0:i]/d[0:i])
    print('主对角项')
    print(d)
    B.ldlfor(a,b)
    for i in range(1,n+1):
        a[i-1,:]=a[i-1,:]/d[i-1]
    B.subbac(a,b)
    print('解向量',b)   
    子程序ldlt
    
    def ldlt(a,d):
      n=a.shape[0]
      for k in range(1,n):
        d[0]=a[0,0]
        if abs(a[k-1,k-1])>1.0e-10:
          for i in range(k+1,n+1):
            x=a[i-1,k-1]/a[k-1,k-1]
            for j in range(k+1,n+1):
              a[i-1,j-1]=a[i-1,j-1]-a[k-1,j-1]*x
            d[i-1]=a[i-1,i-1]
        else:
          print('有0向量在第',k,'行')
    
    子程序ldlfor
    
    def ldlfor(a,b):
      n=a.shape[0]
      for i in range(1,n+1):
        total=b[i-1]
        if i>1:
          for j in range(1,i):
            total=total-a[j-1,i-1]*b[j-1]
        b[i-1]=total/a[i-1,i-1]
    子程序subbac
    
    def subbac(a,b):
      n=a.shape[0]
      for i in range(n,0,-1):
        total=b[i-1]
        if i<n:
          for j in range(i+1,n+1):
            total=total-a[i-1,j-1]*b[j-1]
        b[i-1]=total/a[i-1,i-1]

    程序终端输出结果:
    在这里插入图片描述
    程序结果与计算结果完全一致

    展开全文
  • LDLT分解Python

    千次阅读 2017-11-12 12:35:30
    #!/usr/bin/env python3 ...LDLT ''' import numpy as np def LDLT(amatrix): if len(np.shape(amatrix)) != 2 or np.shape(amatrix)[0]!=np.shape(amatrix)[1]: print("error shape") ret...
  • matlab中矩阵LDLT分解与Cholesky分解矩阵LDLT分解与Cholesky分解:矩阵的LDLT消去函数的程序代码:%矩阵的LDLT分解function [s,l,d]=ldlt(a)s=1;l=0;d=0;%判断矩阵是否对称if a~=a' %矩阵不对称,输出错误信息 s=0;...
  • LU分解、LDLT分解和Cholesky分解

    万次阅读 多人点赞 2014-03-24 13:57:21
    LU分解 概念:假定我们能把矩阵A写成下列两个矩阵相乘的形式:A=LU,其中L为下三角矩阵,U为上三角矩阵。这样我们可以把线性方程组Ax= b写成 Ax= (LU)x = L(Ux) = b。令Ux = y,则原线性方程组Ax = b可首先求解向量y...
  • 对称矩阵的LDLT分解(C语言)

    千次阅读 2012-03-30 18:45:56
    /*对称矩阵的LDLT分解  定理2.2.2:对称矩阵的三角分解:  设A是n阶对称矩阵,若A的各阶顺序主子式均不等于0,则A可以唯一地分解为  A=LDL’ -------------A=LDL’的分解算法------- 参考教材:《数值分析》...
  • 高校计算方法上机作业对矩阵进行LDLT分解及cholesky分解的matlab程序
  • 在线性代数与数值分析中,LU分解是矩阵分解的一种,将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积,有时需要再乘上一个置换矩阵。LU分解可以被视为高斯消元法的矩阵形式。在数值计算上,LU分解经常被用来解...
  • LDLT分解法3. Cholesky分解4. QR分解5.SVD分解6. Jordan 分解 矩阵分解是将矩阵拆解为数个矩阵的乘积,可分为三角分解、满秩分解、QR分解、Jordan分解和SVD(奇异值)分解、非负矩阵分解等,常见的有三种:1)三角...
  • 乔累斯基分解公式简介LLT分解证明具体解法稳定性LDLT分解证明具体解法例子LLT分解LDLT分解引用 矩阵的三角分解是求解线性方程组常用的方法,包括LU分解,LDU分解,杜利特(Doolittle)分解,克劳特(Crout)分解,LLT(乔...
  • 前言 本博客主要介绍在SLAM问题中常常出现的一些线性代数相关的知识,重点是如何采用矩阵分解的方法,求解线性方程组AX=B。主要参考了《计算机...2、LDLT分解与LLT分解(Cholesky分解) 3、QR分解 4、奇异值分...
  • 本博客主要介绍在SLAM问题中常常出现的一些线性代数相关的知识,重点是如何采用矩阵分解的方法,...包括:1、三角分解(LU分解)2、QR分解3、特征值分解4、奇异值分解(SVD分解)5、LDLT分解6、LLT分解(Cholesky分解)
  • holesky分解求dx https://blog.csdn.net/zhouliyang1990/article/details/21952485?t=1501943279852
  • 矩阵分解总结

    千次阅读 2016-12-23 21:27:31
    1. 简介对于线性系统Ax=bAx=b,其解为x=A−1bx=A^{-1}b。但通常情况下,求矩阵的逆十分耗时,应尽量避免。常用的方法包括直接法(Direct Method)和迭代法。   直接法的核心就是矩阵分解(Matrix ...LDLT分解
  • 文章目录L U​ 分解用LU 分解解线性方程组消元法的计算量LU分解的存在性和唯一性对称矩阵的 LDLTL D L^{T}LDLT 分解置换矩阵PA=LUP A=L UPA=LU 分解参考 L U​ 分解 回忆消元法的过程:方阵 AAA 经过初等行变换 ⟶\...
  • 矩阵的LU分解法——Python实现

    万次阅读 多人点赞 2019-04-11 02:31:46
    LU分解、LDLT分解和Cholesky分解 https://blog.csdn.net/zhouliyang1990/article/details/21952485 Doolittle分解法(LU分解)详细分析以及matlab的实现 https://blog.csdn.net/lol_IP/article/details/78491457 ...
  • 超定方程组解法: 最小二乘解法 LDLT分解法解超定方程
  • 矩阵的三角分解是求解线性方程组常用的方法,包括LU分解,LDU分解,杜利特(Doolittle)分解,克劳特(Crout)分解,LLT(乔累斯基Cholesky)分解LDLT(不带平方根乔累斯基)分解等,以及为了满足分解条件又加入行列变换的...
  • 矩阵的三角分解是求解线性方程组常用的方法,包括LU分解,LDU分解,杜利特(Doolittle)分解,克劳特(Crout)分解,LLT(乔累斯基Cholesky)分解LDLT(不带平方根乔累斯基)分解等,以及为了满足分解条件又加入行列变换的...
  • 矩阵的三角分解是求解线性方程组常用的方法,包括LU分解,LDU分解,杜利特(Doolittle)分解,克劳特(Crout)分解,LLT(乔累斯基Cholesky)分解LDLT(不带平方根乔累斯基)分解等,以及为了满足分解条件又加入行列变换的...
  • 并行计算LDLt分解,运用openMP运行速度变慢了,为什么?
  • 矩阵分解

    2020-03-25 13:34:41
    (LDLT 和LL分解合起来称为乔里斯基分解Cholesky decomposition) 奇异值分解 特征分解分解 正交矩阵 正交矩阵:是一个方块矩阵 Q,其元素为实数,而且行向量与列向量皆为正交的单位向量,使得该...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 532
精华内容 212
关键字:

ldlt分解