6.
。
5.Rodrigues
李代数中有三种求导方式：基于指数映射求导、基于BCH公式求导、基于扰动方式求导。
三种求导方式的具体理论及其如何应用并不能一两句话讲清楚，具体可参见相关文献。
这里主要是验证基于指数映射求导和基于BCH公式求导的一致性，OpenCV中Rodrigues是基于指数映射求导(具体可参见源码)，这里实现基于BCH公式求导，具体公式如下图。
以下是详细代码，依赖于C++14、OpenCV4.x和Spdlog。
1 #include <opencv2/opencv.hpp> 2 #include <spdlog/spdlog.h> 3 using namespace std; 4 using namespace cv; 5 6 #ifndef RAD2DEG 7 #define RAD2DEG (180 * 0.3183098861837906715) 8 #define DEG2RAD (3.14159265358979323846 * 0.0055555555555555556) 9 #endif 10 11 Matx33d eulerRot(Matx31d radian) 12 { 13 Matx33d R; 14 double sinR = sin(radian.val[0]); 15 double sinP = sin(radian.val[1]); 16 double sinY = sin(radian.val[2]); 17 double cosR = cos(radian.val[0]); 18 double cosP = cos(radian.val[1]); 19 double cosY = cos(radian.val[2]); 20 21 //RPY indicates: first Yaw aroundZ, second Pitch aroundY, third Roll aroundX 22 R.val[0] = cosY * cosP; R.val[1] = cosY * sinP * sinR - sinY * cosR; R.val[2] = cosY * sinP * cosR + sinY * sinR; 23 R.val[3] = sinY * cosP; R.val[4] = sinY * sinP * sinR + cosY * cosR; R.val[5] = sinY * sinP * cosR - cosY * sinR; 24 R.val[6] = -sinP; R.val[7] = cosP * sinR; R.val[8] = cosP * cosR; 25 return R; 26 } 27 28 static void checkRodrigues(int argc = 0, char** argv = 0) 29 { 30 int N = 999; 31 for (int k = 0; k < N; ++k) 32 { 33 //1.GenerateSimDataAndGT 34 Matx31d degree = Matx31d::randu(-180, 180); 35 Matx33d R = eulerRot(degree * DEG2RAD); 36 Matx31d r; cv::Rodrigues(R, r); 37 Matx31d PW(r.randu(-999, 999)(0), r.randu(-999, 999)(0), r.randu(0, 999)(0)); 38 Matx31d PC = R * PW; 39 40 //2.CalcByExpMap 41 Matx<double, 3, 9> dPCdR; dPCdR << 42 PW(0), PW(1), PW(2), 0, 0, 0, 0, 0, 0, 43 0, 0, 0, PW(0), PW(1), PW(2), 0, 0, 0, 44 0, 0, 0, 0, 0, 0, PW(0), PW(1), PW(2); 45 Mat_<double> dRdr; 46 cv::Rodrigues(r, Matx33d(), dRdr); 47 transpose(dRdr, dRdr); 48 Matx33d dPCdr1 = dPCdR * Matx<double, 9, 3>(dRdr.ptr<double>()); 49 50 //3.CalcByBCH 51 double theta = sqrt(r.val[0] * r.val[0] + r.val[1] * r.val[1] + r.val[2] * r.val[2]); 52 double itheta = 1 / theta; 53 double sn = sin(theta); 54 double cs1 = 1 - cos(theta); 55 Matx31d n(r.val[0] * itheta, r.val[1] * itheta, r.val[2] * itheta); 56 Matx33d nskew(0, -n.val[2], n.val[1], n.val[2], 0, -n.val[0], -n.val[1], n.val[0], 0); 57 Matx33d Jl = itheta * sn * Matx33d::eye() + itheta * cs1 * nskew + (1 - itheta * sn) * n * n.t(); 58 Matx33d skewPC(0, -PC.val[2], PC.val[1], PC.val[2], 0, -PC.val[0], -PC.val[1], PC.val[0], 0); 59 Matx33d dPCdr2 = -skewPC * Jl; 60 61 //4.AnalyzeError 62 double infdPCdr1dPCdr2 = norm(dPCdr1, dPCdr2, NORM_INF); 63 64 //5.PrintError 65 cout << endl << "LoopCount: " << k << endl; 66 if (infdPCdr1dPCdr2 > 1e-9) 67 { 68 cout << endl << "5.1PrintError" << endl; 69 cout << endl << "infdPCdr1dPCdr2: " << infdPCdr1dPCdr2 << endl; 70 if (0) 71 { 72 cout << endl << "5.2PrintDiff" << endl; 73 cout << endl << "dPCdr1: " << endl << dPCdr1 << endl; 74 cout << endl << "dPCdr2: " << endl << dPCdr2 << endl; 75 cout << endl << "5.3PrintOthers" << endl; 76 cout << endl << "PW: " << endl << PW << endl; 77 cout << endl << "PC: " << endl << PC << endl; 78 cout << endl << "r-degree: " << endl << r.t() << endl << degree.t() << endl; 79 } 80 cout << endl << "Press any key to continue" << endl; 81 std::getchar(); 82 } 83 } 84 } 85 86 int main(int argc, char** argv) { checkRodrigues(argc, argv); return 0; }
4.PCAProject
关于PCA的基础知识如下图像所示。
提供的checkPCAProject具有以下目的：验证以上提到的三种计算方法的计算结果与PCAProject的结果的一致性，实现理论验证和接口测试的双重目的。
由于不同的计算方法得到的特征向量可能存在正负号的差别(这对实际应用无影响)，这就导致主分量变换的结果并不一样，所以代码中没有直接判断正向变换的一致性，而是判断重建结果的一致性。
以下是详细代码，依赖于C++14、OpenCV4.x和Spdlog。
1 #include <opencv2/opencv.hpp> 2 #include <spdlog/spdlog.h> 3 using namespace std; 4 using namespace cv; 5 6 #ifndef RandomMM 7 #define RandomMM(min, max) (rand() % ((max) - (min) + 1) + (min)) 8 #endif 9 10 static void checkPCAProject(int argc = 0, char** argv = 0) 11 { 12 int N = 99; 13 for (int k = 0; k < N; ++k) 14 { 15 //1.GenerateSimDataAndGT 16 Mat_<double> X(RandomMM(100, 999), RandomMM(1, 99)); randu(X, -999, 999); 17 Mat_<double> C, O; cv::calcCovarMatrix(X, C, O, COVAR_NORMAL | COVAR_ROWS | COVAR_SCALE); 18 19 //2.CalcByEigenDecompsition 20 Mat_<double> w1, vt1; 21 cv::eigen(C, w1, vt1); 22 23 //3.CalcBySVDCovarMatrix 24 Mat_<double> w2, vt2, u2; 25 cv::SVDecomp(C, w2, u2, vt2, SVD::FULL_UV); 26 27 //4.CalcBySVDCenteredSampleMatrix 28 Mat_<double> w3, vt3, u3; 29 Mat_<double> Xo(X.size()); for (int i = 0; i < X.cols; ++i) Xo.col(i) = X.col(i) - O(i); 30 cv::SVDecomp(Xo, w3, u3, vt3, SVD::FULL_UV); 31 pow(w3, 2, w3); w3 *= (1. / X.rows); 32 33 //5.CalcByPCA 34 Mat_<double> w4, vt4; 35 cv::PCACompute(X, O, vt4, w4, 0); 36 37 //6.PCAProject 38 Mat_<double> Y1 = Xo * vt1.t(); 39 Mat_<double> Y2 = Xo * vt2.t(); 40 Mat_<double> Y3 = Xo * vt3.t(); 41 Mat_<double> Y4 = Xo * vt4.t(); 42 Mat_<double> Y5; cv::PCAProject(X, O, vt4, Y5); 43 44 //7.PACBackProject 45 Mat_<double> Z1 = Y1 * vt1; for (int i = 0; i < X.cols; ++i) Z1.col(i) += O(i); 46 Mat_<double> Z2 = Y2 * vt2; for (int i = 0; i < X.cols; ++i) Z2.col(i) += O(i); 47 Mat_<double> Z3 = Y3 * vt3; for (int i = 0; i < X.cols; ++i) Z3.col(i) += O(i); 48 Mat_<double> Z4 = Y4 * vt4; for (int i = 0; i < X.cols; ++i) Z4.col(i) += O(i); 49 Mat_<double> Z5; cv::PCABackProject(Y5, O, vt4, Z5); 50 51 //8.AnalyzeError 52 double infY1Y5 = norm(Y1, Y5, NORM_INF); 53 double infY2Y5 = norm(Y2, Y5, NORM_INF); 54 double infY3Y5 = norm(Y3, Y5, NORM_INF); 55 double infY4Y5 = norm(Y4, Y5, NORM_INF); 56 double infZ1Z5 = norm(Z1, Z5, NORM_INF); 57 double infZ2Z5 = norm(Z2, Z5, NORM_INF); 58 double infZ3Z5 = norm(Z3, Z5, NORM_INF); 59 double infZ4Z5 = norm(Z4, Z5, NORM_INF); 60 61 //9.PrintError 62 cout << endl << "LoopCount: " << k << endl; 63 if (/*infY1Y5 > 1e-6 || infY2Y5 > 1e-6 || infY3Y5 > 1e-6 || infY4Y5 > 1e-6 || */infZ1Z5 > 1e-6 || infZ2Z5 > 1e-6 || infZ3Z5 > 1e-6 || infZ4Z5 > 1e-6) 64 { 65 cout << endl << "5.1PrintError" << endl; 66 cout << endl << "infY1Y5: " << infY1Y5 << endl; 67 cout << endl << "infY2Y5: " << infY2Y5 << endl; 68 cout << endl << "infY3Y5: " << infY3Y5 << endl; 69 cout << endl << "infY4Y5: " << infY4Y5 << endl; 70 cout << endl << "infZ1Z5: " << infZ1Z5 << endl; 71 cout << endl << "infZ2Z5: " << infZ2Z5 << endl; 72 cout << endl << "infZ3Z5: " << infZ3Z5 << endl; 73 cout << endl << "infZ4Z5: " << infZ4Z5 << endl; 74 if (0) 75 { 76 cout << endl << "5.2PrintDiff" << endl; 77 cout << endl << "Y1: " << endl << Y1 << endl; 78 cout << endl << "Y2: " << endl << Y2 << endl; 79 cout << endl << "Y3: " << endl << Y3 << endl; 80 cout << endl << "Y4: " << endl << Y4 << endl; 81 cout << endl << "Y5: " << endl << Y5 << endl; 82 cout << endl; 83 cout << endl << "Z1: " << endl << Z1 << endl; 84 cout << endl << "Z2: " << endl << Z2 << endl; 85 cout << endl << "Z3: " << endl << Z3 << endl; 86 cout << endl << "Z4: " << endl << Z4 << endl; 87 cout << endl << "Z5: " << endl << Z5 << endl; 88 cout << endl << "5.3PrintOthers" << endl; 89 cout << endl << "C: " << endl << C << endl; 90 cout << endl << "X: " << endl << X << endl; 91 } 92 cout << endl << "Press any key to continue" << endl; 93 std::getchar(); 94 } 95 } 96 } 97 98 int main(int argc, char** argv) { checkPCAProject(argc, argv); return 0; }
3.eigen
仅方阵且可相似对角化才能进行特征值分解，关于相似对角化有以下：
(1)充要条件：有n个线性无关的特征向量
(2)充分条件：k重特征值有k个线性无关特征向量(这样就肯定有n个线性无关的特征向量)
(3)充分条件：有n个不同特征值(因为有定理表述不同特征值对应的特征向量线性无关)
(4)充分条件：为正规矩阵(因为正规矩阵要么有n个不同的特征值要么k重特征值对应k个线性无关的特征向量)
提供的checkEigen具有以下目的：验证特征值分解eigen、奇异值分解SVD、主成分分析PCA三者结果的一致性，或者说基于后两者验证前者的正确性。
以下是详细代码，依赖于C++14、OpenCV4.x和Spdlog。
1 #include <opencv2/opencv.hpp> 2 #include <spdlog/spdlog.h> 3 using namespace std; 4 using namespace cv; 5 6 #ifndef RandomMM 7 #define RandomMM(min, max) (rand() % ((max) - (min) + 1) + (min)) 8 #endif 9 10 static void checkEigen(int argc = 0, char** argv = 0) 11 { 12 int N = 99; 13 for (int k = 0; k < N; ++k) 14 { 15 //1.GenerateSimDataAndGT 16 Mat_<double> X(RandomMM(111, 999), RandomMM(11, 99)); cv::randu(X, -999, 999); 17 Mat_<double> C, O; cv::calcCovarMatrix(X, C, O, COVAR_NORMAL | COVAR_ROWS | COVAR_SCALE); 18 19 //2.CalcByEigenDecompsition 20 Mat_<double> w1, vt1; 21 cv::eigen(C, w1, vt1); 22 23 //3.CalcBySVDCovarMatrix 24 Mat_<double> w2, vt2, u2; 25 cv::SVDecomp(C, w2, u2, vt2, SVD::FULL_UV); 26 27 //4.CalcBySVDCenteredSampleMatrix 28 Mat_<double> w3, vt3, u3; 29 Mat_<double> Xo(X.size()); for (int i = 0; i < X.cols; ++i) Xo.col(i) = X.col(i) - O(i); 30 cv::SVDecomp(Xo, w3, u3, vt3, SVD::FULL_UV); 31 pow(w3, 2, w3); w3 *= (1. / X.rows); 32 33 //5.CalcByPCA 34 Mat_<double> w4, vt4; 35 cv::PCACompute(X, O, vt4, w4); 36 37 //6.AnalyzeError 38 double infvt1vt2 = norm(cv::abs(vt1), cv::abs(vt2), NORM_INF); 39 double infvt1vt3 = norm(cv::abs(vt1), cv::abs(vt3), NORM_INF); 40 double infvt1vt4 = norm(cv::abs(vt1), cv::abs(vt4), NORM_INF); 41 double infw1w2 = norm(w1, w2, NORM_INF); 42 double infw1w3 = norm(w1, w3, NORM_INF); 43 double infw1w4 = norm(w1, w4, NORM_INF); 44 45 //7.PrintError 46 cout << endl << "LoopCount: " << k << endl; 47 if (infvt1vt2 > 1e-6 || infvt1vt3 > 1e-6 || infvt1vt4 > 1e-6 || infw1w2 > 1e-6 || infw1w3 > 1e-6 || infw1w4 > 1e-6) 48 { 49 cout << endl << "5.1PrintError" << endl; 50 cout << endl << "infvt1vt2: " << infvt1vt2 << endl; 51 cout << endl << "infvt1vt3: " << infvt1vt3 << endl; 52 cout << endl << "infvt1vt4: " << infvt1vt4 << endl; 53 cout << endl << "infw1w2: " << infw1w2 << endl; 54 cout << endl << "infw1w3: " << infw1w2 << endl; 55 cout << endl << "infw1w4: " << infw1w2 << endl; 56 if (0) 57 { 58 cout << endl << "5.2PrintDiff" << endl; 59 cout << endl << "vt1: " << endl << vt1 << endl; 60 cout << endl << "vt2: " << endl << vt2 << endl; 61 cout << endl << "vt3: " << endl << vt3 << endl; 62 cout << endl << "vt4: " << endl << vt4 << endl; 63 cout << endl << "w1: " << endl << w1 << endl; 64 cout << endl << "w2: " << endl << w2 << endl; 65 cout << endl << "w3: " << endl << w3 << endl; 66 cout << endl << "w4: " << endl << w4 << endl; 67 cout << endl << "5.3PrintOthers" << endl; 68 cout << endl << "C: " << endl << C << endl; 69 cout << endl << "X: " << endl << X << endl; 70 } 71 cout << endl << "Press any key to continue" << endl; 72 std::getchar(); 73 } 74 } 75 } 76 77 int main(int argc, char** argv) { checkEigen(argc, argv); return 0; }
2.calcCovarMatrix
关于期望和协方差的基础知识如下图像所示。
提供的checkCalcCovarMatrix具有以下目的：验证公式计算的结果与calcCovarMatrix的结果的一致性，实现理论验证和接口测试的双重目的。
以下是详细代码，依赖于C++14、OpenCV4.x和Spdlog。
1 #include <opencv2/opencv.hpp> 2 #include <spdlog/spdlog.h> 3 using namespace std; 4 using namespace cv; 5 6 #ifndef RandomMM 7 #define RandomMM(min, max) (rand() % ((max) - (min) + 1) + (min)) 8 #endif 9 10 static void checkCalcCovarMatrix(int argc = 0, char** argv = 0) 11 { 12 int N = 999; 13 for (int k = 0; k < N; ++k) 14 { 15 //1.GenerateSimDataAndGT 16 Mat_<double> X(RandomMM(111,999), RandomMM(111,999)); randu(X, -999, 999); 17 18 //2.CalcByOpenCV 19 Mat_<double> C1, O1; 20 cv::calcCovarMatrix(X, C1, O1, COVAR_NORMAL | COVAR_ROWS | COVAR_SCALE); 21 22 //3.CalcByTheory 23 Mat_<double> O2(1, X.cols); for (int j = 0; j < X.cols; ++j) O2(j) = mean(X.col(j))[0]; 24 Mat_<double> C2; mulTransposed(X, C2, true, O2, 1. / X.rows); 25 //Mat_<double> C3 = 1. / X.rows * X.t() * X - O2.t() * O2; C2 = C3; 26 27 //4.AnalyzeError 28 double infO1O2 = norm(O1, O2, NORM_INF); 29 double infC1C2 = norm(C1, C2, NORM_INF); 30 31 //5.PrintError 32 cout << endl << "LoopCount: " << k << endl; 33 if (infO1O2 > 0 || infC1C2 > 0) 34 { 35 cout << endl << "5.1PrintError" << endl; 36 cout << endl << "infO1O2: " << infO1O2 << endl; 37 cout << endl << "infC1C2: " << infC1C2 << endl; 38 if (0) 39 { 40 cout << endl << "5.2PrintDiff" << endl; 41 cout << endl << "O1: " << endl << O1 << endl; 42 cout << endl << "O2: " << endl << O2 << endl; 43 cout << endl << "C1: " << endl << C1 << endl; 44 cout << endl << "C2: " << endl << C2 << endl; 45 cout << endl << "5.3PrintOthers" << endl; 46 cout << endl << "X: " << endl << X << endl; 47 } 48 cout << endl << "Press any key to continue" << endl; 49 std::getchar(); 50 } 51 } 52 } 53 54 int main(int argc, char** argv) { checkCalcCovarMatrix(argc, argv); return 0; }
1.SVDecomp
任意矩阵A都能进行奇异值分解，且不论A的是否为方阵或是否满秩(包括行满秩、列满秩及全满秩)，A的SVD形式都可以归结为两种：完整型和缩减型。
假设Amn是m行n列的矩阵且秩为r，则：
(1)完整型：Amn=Umm*Wmn*Trans(Vnn)
(2)缩减型：Amn=Umr*Wrr*Trans(Vnr)
根据A是否满秩，可根据以上两种形式衍生出多种形式，但无非都是r或取m或n等变化，典型地，满秩方阵的两种形式一致。
SVDecomp正是提供了以上两种分解形式，通过设置flags可得到不同的分解形式，flags取值如下：
(1)NO_UV：仅返回奇异值向量W，返回空U和空V。
(2)FULL_UV：按完整型进行SVD，对满秩方阵无影响，对欠秩矩阵和非方阵有作用。
(3)MODIFY_A：暂无作用。
默认是缩减型SVD。需要注意的是无论flags取什么值，返回的W都是维度相同的向量且维度等于行数和列数的较小者，只是对于欠秩矩阵，W最后几维的数值为0。
提供的checkSVDecomp具有以下目的：
(1)如何使用SVDecomp。
(2)测试SVDecomp的正确性：通过对比分析原始矩阵和重建矩阵来测试。
以下是详细代码，依赖于C++14、OpenCV4.x和Spdlog。
1 #include <opencv2/opencv.hpp> 2 #include <spdlog/spdlog.h> 3 using namespace std; 4 using namespace cv; 5 6 static void checkSVDecomp(int argc = 0, char** argv = 0) 7 { 8 int N = 999; 9 for (int k = 0; k < N / 2; ++k) 10 { 11 //1.GenerateSimDataAndGT 12 Mat_<double> A(5, 3); cv::randu(A, -999, 999); 13 14 //2.DecomposeByShortSVD 15 Mat_<double> U1, W1, VT1; 16 cv::SVDecomp(A, W1, U1, VT1); 17 18 //3.DecomposeByFullSVD 19 Mat_<double> U2, W2, VT2; 20 cv::SVDecomp(A, W2, U2, VT2, SVD::FULL_UV); 21 22 //4.AnalyzeError 23 double infA0A1 = norm(U1 * Matx33d::diag(W1) * VT1, A, NORM_INF); 24 Mat_<double> W2_(A.size(), 0.); for (int k = 0; k < 3; ++k) W2_(k, k) = W2(k); 25 double infA0A2 = norm(U2 * W2_ * VT2, A, NORM_INF); 26 27 //5.PrintError 28 cout << endl << "LoopCount: " << k << endl; 29 if (infA0A1 > 1E-6 || infA0A2 > 1E-6) 30 { 31 cout << endl << "5.1PrintError" << endl; 32 cout << endl << "infA0A1: " << infA0A1 << endl; 33 cout << endl << "infA0A2: " << infA0A2 << endl; 34 cout << endl << "5.2PrintDiff" << endl; 35 cout << endl << "U1:" << endl << U1 << endl; 36 cout << endl << "U2:" << endl << U2 << endl; 37 cout << endl << "W1:" << endl << W1 << endl; 38 cout << endl << "W2:" << endl << W2 << endl; 39 cout << endl << "VT1:" << endl << VT1 << endl; 40 cout << endl << "VT2:" << endl << VT2 << endl; 41 cout << endl << "5.3PrintOthers" << endl; 42 cout << endl << "Press any key to continue" << endl; 43 std::getchar(); 44 } 45 } 46 for (int k = N / 2; k < N; ++k) 47 { 48 //1.GenerateSimDataAndGT 49 Mat_<double> A(3, 5); cv::randu(A, -999, 999); 50 51 //2.DecomposeByShortSVD 52 Mat_<double> U1, W1, VT1; 53 cv::SVDecomp(A, W1, U1, VT1); 54 55 //3.DecomposeByFullSVD 56 Mat_<double> U2, W2, VT2; 57 cv::SVDecomp(A, W2, U2, VT2, SVD::FULL_UV); 58 59 //4.AnalyzeError 60 double infA0A1 = norm(U1 * Matx33d::diag(W1) * VT1, A, NORM_INF); 61 Mat_<double> W2_(A.size(), 0.); for (int k = 0; k < 3; ++k) W2_(k, k) = W2(k); 62 double infA0A2 = norm(U2 * W2_ * VT2, A, NORM_INF); 63 64 //5.PrintError 65 cout << endl << "LoopCount: " << k << endl; 66 if (infA0A1 > 1E-6 || infA0A2 > 1E-6) 67 { 68 cout << endl << "5.1PrintError" << endl; 69 cout << endl << "infA0A1: " << infA0A1 << endl; 70 cout << endl << "infA0A2: " << infA0A2 << endl; 71 cout << endl << "5.2PrintDiff" << endl; 72 cout << endl << "U1:" << endl << U1 << endl; 73 cout << endl << "U2:" << endl << U2 << endl; 74 cout << endl << "W1:" << endl << W1 << endl; 75 cout << endl << "W2:" << endl << W2 << endl; 76 cout << endl << "VT1:" << endl << VT1 << endl; 77 cout << endl << "VT2:" << endl << VT2 << endl; 78 cout << endl << "5.3PrintOthers" << endl; 79 cout << endl << "Press any key to continue" << endl; 80 std::getchar(); 81 } 82 } 83 } 84 85 int main(int argc, char** argv) { checkSVDecomp(argc, argv); return 0; }