矩阵迹运算介绍及C++/OpenCV/Eigen的三种实现

矩阵迹运算返回的是矩阵对角元素的和：

迹运算因为很多原因而有用。若不使用求和符号，有些矩阵运算很难描述，而通过矩阵乘法和迹运算符号，可以清楚地表示。例如，迹运算提供了另一种描述矩阵Frobenius范数的方式：

用迹运算表示表达式，我们可以使用很多有用的等式巧妙地处理表达式。例如，迹运算在转置运算下是不变的：Tr(A)=Tr(A^T).

多个矩阵相乘得到的方阵的迹，和将这些矩阵中的最后一个挪到最前面之后相乘的迹是相同的。当然，我们需要考虑挪到之后矩阵乘积依然定义良好：Tr(ABC)=Tr(CAB)=Tr(BCA).

或者更一般地，

即使循环置换后矩阵乘积得到的矩阵形状变了，迹运算的结果依然不变。例如，假设矩阵A∈R^m*n,矩阵B∈R^n*m，我们可以得到Tr(AB)=Tr(BA)，尽管AB∈R^m*m和BA∈R^n*n。

另一个有用的事实是标量在迹运算后仍然是它自己：a=Tr(a)。

在线性代数中，一个n*n的矩阵A的迹(或迹数)，是指A的主对角线(从左上方至右下方的对角线)上各个元素的总和，一般记作tr(A)或Sp(A)：tr(A)=A_1,1+A_2,2+…+A_n,n

其中A_i,j代表矩阵的第i行第j列上的元素的值。一个矩阵的迹是其特征值得总和(按代数重数计算)。迹的英文为trace。

迹运算性质：

(1)、对于任两个n*n的矩阵A、B和标量r，都有：tr(A+B)=tr(A)+tr(B); tr(r·A)=r·tr(A)；

(2)、由于一个矩阵A的转置矩阵A^T的主对角线元素和原来矩阵的主对角元素是一样的，所以任意一个矩阵和其转置矩阵都会有相同的迹：tr(A)=tr(A^T)；

(3)、矩阵乘积的迹数：设A是一个n*m矩阵，B是个m*n矩阵，则tr(AB)=tr(BA)，其中AB是一个n*n矩阵，而BA是一个m*m矩阵；

(4)、如果A和B都是n*n的方形矩阵，那么它们的乘积AB和BA也会是方形矩阵。因此，利用这个结果，可以推导出：计算若干个同样大小的方形矩阵的乘积的迹数时，可以循环改变乘积中方形矩阵相乘的顺序，而最终的结果不变。例如，有三个方形矩阵A、B和C，则tr(ABC)=tr(BCA)=tr(CAB)，但是要注意：tr(ABC)≠tr(ACB)；

    更一般地，乘积中的矩阵不一定要是方形矩阵，只要某一个循环改变后的乘积依然存在，那么得到的迹数依然会和原来的迹数相同。另外，如果A、B和C是同样大小的方阵而且还是对称矩阵的话，那么其乘积的迹数不只在循环置换下不会改变，而且在所有的置换下都不会改变：tr(ABC)=tr(BCA)=tr(CAB)=tr(ACB)=tr(CBA)=tr(BAC)；

    (5)、迹数的相似不变性：迹数拥有相似不变性。如果矩阵A和B相似的话，它们会有相同的迹；矩阵A和B相似也就是说存在可逆矩阵P,使得B=PAP^-1.

以上内容摘自 《深度学习中文版》 和 维基百科

以下是分别采用C++和OpenCV实现的求迹运算code：

#include "funset.hpp"
#include <math.h>
#include <iostream>
#include <string>
#include <vector>
#include <opencv2/opencv.hpp>
#include "common.hpp"

// ================================= 求矩阵的迹 =================================
template<typename _Tp>
_Tp trace(const std::vector<std::vector<_Tp>>& mat)
{
	_Tp ret{ (_Tp)0 };
	int nm = std::min(mat.size(), mat[0].size());

	for (int i = 0; i < nm; ++i) {
		ret += mat[i][i];
	}

	return ret;
}

int test_trace()
{
	std::vector<std::vector<float>> vec{ { 1.2f, 2.5f, 5.6f, -2.5f },
					{ -3.6f, 9.2f, 0.5f, 7.2f },
					{ 4.3f, 1.3f, 9.4f, -3.4f } };
	const int rows{ 3 }, cols{ 4 };

	fprintf(stderr, "source matrix:\n");
	print_matrix(vec);

	float tr = trace(vec);
	fprintf(stderr, "\nc++ implement trace: %f\n", tr);

	cv::Mat mat(rows, cols, CV_32FC1);
	for (int y = 0; y < rows; ++y) {
		for (int x = 0; x < cols; ++x) {
			mat.at<float>(y, x) = vec.at(y).at(x);
		}
	}

	cv::Scalar scalar = cv::trace(mat);
	fprintf(stderr, "\nopencv implement trace: %f\n", scalar.val[0]);

	return 0;
}

执行结果如下：

以下是采用Eigen实现的迹运算code：

#include "funset.hpp"
#include <math.h>
#include <iostream>
#include <vector>
#include <string>
#include <opencv2/opencv.hpp>
#include <Eigen/Dense>
#include "common.hpp"

int test_trace()
{
	std::vector<std::vector<float>> vec{ { 1.2f, 2.5f, 5.6f, -2.5f },
					{ -3.6f, 9.2f, 0.5f, 7.2f },
					{ 4.3f, 1.3f, 9.4f, -3.4f } };
	const int rows{ 3 }, cols{ 4 };

	std::vector<float> vec_;
	for (int i = 0; i < rows; ++i) {
		vec_.insert(vec_.begin() + i * cols, vec[i].begin(), vec[i].end());
	}
	Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> m(vec_.data(), rows, cols);

	fprintf(stderr, "source matrix:\n");
	std::cout << m << std::endl;

	float tr = m.trace();
	fprintf(stderr, "\nEigen implement trace: %f\n", tr);

	return 0;
}

执行结果如下：

由以上结果可见：C++、OpenCV、Eigen实现结果是一致的。

GitHub：

https://github.com/fengbingchun/NN_Test

https://github.com/fengbingchun/Eigen_Test