Eigen学习日志3

#include "iostream"
#include <Eigen/Dense>
#include <Eigen/Core>
#include <Eigen/Eigen>

using namespace Eigen;

//关于slam算法为什么要进行矩阵分解 https://blog.csdn.net/wangshuailpp/article/details/80209863?utm_medium=distribute.pc_relevant.none-task-blog-2~default~baidujs_baidulandingword~default-4.no_search_link&spm=1001.2101.3001.4242.3


//向量范数，矩阵范数
void NormComputations();
//基础线性方程的求解AX=b A为矩阵，b为向量，x为待求解向量
void BasicLinearSolving();
//矩阵特征值，特征向量计算
void ComputingEigenvaluesAndEigenvectors();

void NormComputations(){

    VectorXf v(2);
    MatrixXf m(2,2);

    v << -3,-2;
    m << 1,-2,-3,4;
    std::cout<<"v:"<<std::endl<<v<<std::endl;
    std::cout<<"m:"<<std::endl<<m<<std::endl;
    //向量范数
    //squaredNorm()，未进行开方的l2范数即 ||x||2^2 = x1^2+x2^2+.......
    std::cout << "v.squaredNorm() = " << v.squaredNorm() << std::endl;
    //norm(),返回的是l2范数
    std::cout << "v.norm() = " << v.norm() << std::endl;
    //lpNorm<1>()，返回的是l1范数，即向量各元素绝对值之和
    std::cout << "v.lpNorm<1>() = " << v.lpNorm<1>() << std::endl;
    //lpNorm<Infinity>()，返回各元素中绝对值最大值
    std::cout << "v.lpNorm<Infinity>() = " << v.lpNorm<Infinity>() << std::endl;
    std::cout << std::endl;
    //矩阵范数
    //squaredNorm()，未进行开方的l2范数即 ||x||2^2 = x1^2+x2^2+.......
    std::cout << "m.squaredNorm() = " << m.squaredNorm() << std::endl;
    //norm(),返回的是l2范数
    std::cout << "m.norm() = " << m.norm() << std::endl;
    //lpNorm<1>()，返回的是l1范数，即向量各元素绝对值之和
    std::cout << "m.lpNorm<1>() = " << m.lpNorm<1>() << std::endl;
    //lpNorm<Infinity>()，返回各元素中绝对值最大值
    std::cout << "m.lpNorm<Infinity>() = " << m.lpNorm<Infinity>() << std::endl;

    //矩阵的1-范数1-norm：列和范数，即所有矩阵列向量绝对值之和的最大值
    //计算方式如下
    std::cout<<"1-norm(m)="<<m.cwiseAbs().colwise().sum().maxCoeff()<< " == " << m.colwise().lpNorm<1>().maxCoeff() << std::endl;
    //矩阵的无穷范数：行和范数，即所有矩阵行向量绝对值之和的最大值
    //计算方式如下
    std::cout << "infty-norm(m) = " << m.cwiseAbs().rowwise().sum().maxCoeff()<< " == " << m.rowwise().lpNorm<1>().maxCoeff() << std::endl;
    /*
     *  v:
        -3
        -2
        m:
         1 -2
        -3  4
        v.squaredNorm() = 13
        v.norm() = 3.60555
        v.lpNorm<1>() = 5
        v.lpNorm<Infinity>() = 3

        m.squaredNorm() = 30
        m.norm() = 5.47723
        m.lpNorm<1>() = 10
        m.lpNorm<Infinity>() = 4
        1-norm(m)=6 == 6
        infty-norm(m) = 7 == 7
     */


}
void BasicLinearSolving(){
    //线性方程求解基础知识

    //对于Ax = 0而言
    /*
     * 对于A为n*n的矩阵而言
     * 若|A|不为0，即A的秩为n时r(A)=n，该方程有唯一解
     * 若|A|= 0，则r(A)<n,该方程有无数解
     *
     * 对于A为m*n的矩阵而言
     * 若r(A)<n,则该方程有无数解
     * 若r(A)=n,则该方程有唯一解
     */

    //对于Ax = b而言
    /*
     * 对于A为m*n的矩阵而言
     * 1.该方程无解的充要条件为 r(A) < r(A,b)
     * 2.该方程有唯一解的充要条件为 r(A) = r(A,b) = n
     * 3.该方程有无数解的充要条件为 r(A) = r(A,b) < n
     */

    //正交矩阵
    /*
     * 如果A*A^T=E（E为单位矩阵，A^T表示“矩阵A的转置矩阵”)或A^T*A=E，则n阶实矩阵A称为正交矩阵
     * 正交矩阵性质
     * 若B为正交矩阵、则B必定为可逆矩阵
     * 则B^(-1) = B^(T)
     *
     */

    //相似矩阵：相似矩阵是指存在相似关系的矩阵。设A，B为n阶矩阵，如果有n阶可逆矩阵P存在，
    // 使得P^(-1)AP=B
    // 则称矩阵A和B相似
    /*
     *  相似矩阵性质
     *  对于
        设A，B和C是任意同阶方阵，则有：
        （1）反身性：A~ A
        （2）对称性：若A~ B，则 B~ A
        （3）传递性：若A~ B，B~ C，则A~ C
        （4）若A~ B，则r(A)=r(B)，|A|=|B|，tr（A）=tr（B）。
        （5）若A~ B，且A可逆，则B也可逆，且B~ A。
        （6）若A~ B，则A与B
            两者的秩相等；
            两者的行列式值相等；
            两者的迹数相等；
            两者拥有同样的特征值，尽管相应的特征向量一般不同；
            两者拥有同样的特征多项式；
            两者拥有同样的初等因子。
        （7）若A与对角矩阵相似，则称A为可对角化矩阵，若n阶方阵A有n个线性无关的特征向量，则称A为单纯矩阵。
        （8）相似矩阵具有相同的可逆性，当它们可逆时，则它们的逆矩阵也相似
     */


    //实对称矩阵：n阶矩阵A，其矩阵的元素都为实数，且矩阵A的转置等于其本身，则A为实对称矩阵
    /*
     * 相似对角化：如果一个方阵A相似于对角矩阵，也就是说存在一个可逆矩阵P使得P^(-1)AP是对角矩阵，则A就被称为可以相似对角化的
     *
     * 实对称矩阵性质
     * 1.实对称矩阵A的不同特征值对应的特征向量是正交的。
       2.实对称矩阵A的特征值都是实数。
       3.n阶实对称矩阵A必可相似对角化，且相似对角阵上的元素即为矩阵本身特征值。
         即存在可逆矩阵P 使得P^(-1)AP=M（M为对角矩阵），则M的对角线元素即为A矩阵本身特征值
         且 A = PMP^(-1)
       4.若A具有k重特征值λ　必有k个线性无关的特征向量，或者说秩r(λE-A)必为n-k，其中E为单位矩阵。
       5.实对称矩阵A一定可正交相似对角化。
     */

    //直接求逆的方法求解线性方程组
    Matrix<float,3,3> m_33f;
    m_33f<<3,1,1,
           1,3,1,
           1,1,3;
    std::cout<<"m_33f:"<<std::endl<<m_33f<<std::endl;
    Vector3f b_3f;
    b_3f<<0,3,2;
    std::cout<<"b_3f:"<<std::endl<<b_3f<<std::endl;

    Matrix<float,3,3> m_inverse_33f = m_33f.inverse();
    std::cout<<"m_inverse_33f:"<<std::endl<<m_inverse_33f<<std::endl;
    Vector3f x_m33f = m_inverse_33f*b_3f;
    std::cout<<" x_m33f = "<<std::endl<<x_m33f<<std::endl;

    /*
     *  m_33f:
        3 1 1
        1 3 1
        1 1 3
        b_3f:
        0
        3
        2
        m_inverse_33f:
         0.4 -0.1 -0.1
        -0.1  0.4 -0.1
        -0.1 -0.1  0.4
         x_m33f =
        -0.5
           1
         0.5
     */

    //QR分解 即对于n阶可逆矩阵（非奇异矩阵|A|不为0）A（nxn），可以分解成正交矩阵Q（nxn）和非奇异（可逆）上三角矩阵R（nxn）的乘积
    //对于A（mxn）矩阵，若r(A)=n,则可以分解成正交矩阵Q（mxn）和非奇异（可逆）上三角矩阵R（nxn）的乘积
    //即A = QR，则称该分解为QR分解 详解见 https://blog.csdn.net/lhxez6868/article/details/100165447
    //QR分解求解步骤 https://blog.csdn.net/wangshuailpp/article/details/80209863?utm_medium=distribute.pc_relevant.none-task-blog-2~default~baidujs_baidulandingword~default-4.no_search_link&spm=1001.2101.3001.4242.3

    //使用QR分解求解 matrix_33f*x=v_3f
    Matrix<float,3,3> matrix_33f;
    matrix_33f<<3,1,1,
                1,3,1,
                1,1,3;
    std::cout<<"matrix_33f:"<<std::endl<<matrix_33f<<std::endl;
    Vector3f v_3f;
    v_3f<<0,3,2;
    std::cout<<"v_3f:"<<std::endl<<v_3f<<std::endl;

    Vector3f x = matrix_33f.colPivHouseholderQr().solve(v_3f);
    std::cout<<" x = "<<std::endl<<x<<std::endl;

    //在此示例中，colPivHouseholderQr（）方法返回类ColPivHouseholderQR的对象。
    //由于此处的矩阵类型为Matrix3f，因此该行可能已替换为：
    //ColPivHouseholderQR <Matrix3f> dec（matrix_33f）;
    //Vector3f x = dec.solve（v_3f）;
    //在这里，ColPivHouseholderQR是具有选择列主元功能的QR分解。
    //这是本教程的一个不错的折衷方案，因为它适用于所有矩阵，而且速度非常快。
    //这是一些其他分解表，您可以根据矩阵和要进行的权衡选择：


    // Decomposition	                         Method	                 Requirements on the matrix	             Speed(small-to-medium)	                Speed (large)	                Accuracy
    // PartialPivLU	                         partialPivLu()	                    Invertible	                            ++	                                ++	                           +
    // FullPivLU	                           fullPivLu()	                        None	                            -	                               - -	                          +++
    // HouseholderQR	                     householderQr()	                    None	                            ++	                                ++                             +
    // ColPivHouseholderQR	              colPivHouseholderQr()	                    None	                            +	                                -	                          +++
    // FullPivHouseholderQR	              fullPivHouseholderQr()	                None                              	-	                                - -	                          +++
    // CompleteOrthogonalDecomposition	completeOrthogonalDecomposition()	        None	                            +	                                -	                          +++
    // LLT	                                      llt()	                        Positive definite	                    +++	                                +++                            +
    // LDLT	                                      ldlt()	             Positive or negative semidefinite	            +++	                                +	                          ++
    // BDCSVD	                                  bdcSvd()	                        None	                            -	                                -	                          +++
    // JacobiSVD	                              jacobiSvd()	                    None	                            -	                                - - -	                      +++
    // 当矩阵是正定矩阵时可以采用ldlt（）方法求解
    //所有这些分解都提供了一个solve（）方法，该方法与上述示例一样。用来求解

    //确认解是否有效，通过计算误差来判断
    double relative_error = (matrix_33f*x-v_3f).norm()/v_3f.norm();
    std::cout<<" relative_error = "<<std::endl<<relative_error<<std::endl;

    //另一种求解误差方式,
    //bool a_solution_exists = （A*result）.isApprox（b， precision）;
    //如果解有效则 isApprox（）返回true，否则为false
    bool solver_exit = v_3f.isApprox(matrix_33f*x);
    std::cout<<solver_exit<<std::endl;
    /*
     *  matrix_33f:
        3 1 1
        1 3 1
        1 1 3
        v_3f:
        0
        3
        2
         x =
        -0.5
           1
         0.5
         relative_error =
        7.39305e-08
        1
     */

}

void ComputingEigenvaluesAndEigenvectors(){
     //特征值，特征向量，方阵的秩
     /*
      * 设: 方阵阶数为n, 特征值个数为k, 其中有i重特征值 λi, 单个λ对应的无关特征向量个数为t, 方阵的秩为r
      * 1.特征值
      * a.特征值个数k（包括重根） 与 方阵的阶数n相等;
      * b.特征值个数k>=所有无关特征向量数之和 (因为i重特征值λi最多有i个线性无关的特征向量);
      * c.特征值个数k与方阵的秩无关。
      *
      * 2. 特征向量
      * a.单个λ对应的无关特征向量个数t与方阵的秩r没有什么直接的关系，它们都小于等于方阵阶数n
      *
      * 3. 方阵的秩
      * a.方阵的秩r与它的特征值λi=0 的重数i有关
      *   i.当方阵A可以相似对角化时（这里自然是包括了方阵A为实对称矩阵的情况），i=n-r.
      *     因为A ~^（对角矩阵）,  所以r(A)=r(^)。 此时若r(A)=r(^)=r, 意味着对角阵有r个不为零的特征值，
      *     即A也有r个不为零的特征值，进而得到A 有n-r 重特征值：λi=0
      *   ii.当方阵A不可相似对角化时，i>=n-r.
      *      首先，对于i重特征值λi最多有i个线性无关的特征向量，
      *      反过来说，同一特征值λi对应的线性无关的特征向量个数（设为s）<= i.
      *      对于λi=0，有r(0E-A)=r(-A)=r(A)=r,所以，λi的线性无关特征向量个数s=n-r, 根据上一行的说法就有，λi=0 的重数 s<=i, 即 i>= n-r.
      */
     //template<typename _MatrixType>
     //class Eigen::SelfAdjointEigenSolver< _MatrixType >  该类用于计算 实对称矩阵 的特征值和特征向量

     //matrix_33f是实对称矩阵
     Matrix3f matrix_33f;
     matrix_33f<<1,3,2,
                 3,1,5,
                 2,5,1;
     std::cout<<" matrix_33f = "<<std::endl<<matrix_33f<<std::endl;
     SelfAdjointEigenSolver<Matrix3f> eigensolver(matrix_33f);
     if (eigensolver.info()!=Eigen::Success)
         abort();
     std::cout << "The eigenvalues of matrix_33f are:\n"
          << eigensolver.eigenvalues() << std::endl;
     std::cout << "Here's a matrix whose columns are eigenvectors of matrix_33f \n"
          << "corresponding to these eigenvalues:\n"
          << eigensolver.eigenvectors() << std::endl;

     //template<typename _MatrixType>
     //class Eigen::EigenSolver< _MatrixType > 该类计算一般矩阵的特征值和特征向量。
     Matrix<double,3,3> matrix_33d;
     matrix_33d<<1,3,2,
                 6,1,4,
                 3,5,2;
     std::cout<<" matrix_33d = "<<std::endl<<matrix_33d<<std::endl;
     EigenSolver<Matrix3d> solver(matrix_33d);
     if (solver.info()!=Eigen::Success)
         abort();
     std::cout << "The eigenvalues of matrix_33d are:\n"
               << solver.eigenvalues() << std::endl;
     std::cout << "Here's a matrix whose columns are eigenvectors of matrix_33d \n"
               << "corresponding to these eigenvalues:\n"
               << solver.eigenvectors() << std::endl;
     /*
      * matrix_33f =
        1 3 2
        3 1 5
        2 5 1
        The eigenvalues of matrix_33f are:
         -4.12839
        -0.710691
          7.83909
        Here's a matrix whose columns are eigenvectors of matrix_33f
        corresponding to these eigenvalues:
         -0.17831 -0.869143  0.461298
         0.738394  0.191674  0.646557
        -0.650369  0.455907  0.607593
         matrix_33d =
        1 3 2
        6 1 4
        3 5 2
        The eigenvalues of matrix_33d are:
         (8.86025,0)
        (-1.07277,0)
        (-3.78749,0)
        Here's a matrix whose columns are eigenvectors of matrix_33d
        corresponding to these eigenvalues:
        (-0.409941,0)  (0.484108,0) (-0.269099,0)
        (-0.642422,0)  (0.228599,0)  (0.793414,0)
        (-0.647489,0)  (-0.84462,0) (-0.545966,0)

        进程已结束，退出代码为 0

      */


}



int main(){

    NormComputations();
    BasicLinearSolving();
    ComputingEigenvaluesAndEigenvectors();
    return  0;

}