伴随矩阵介绍及C++实现

在线性代数中，一个方形矩阵的伴随矩阵是一个类似于逆矩阵的概念。如果矩阵可逆，那么它的逆矩阵和它的伴随矩阵之间只差一个系数。然而，伴随矩阵对不可逆的矩阵也有定义，并且不需要用到除法。

设R是一个交换环(在抽象代数之分支环论中，一个交换环(commutative ring)是乘法运算满足交换律的环)，A是一个以R中元素为系数的n*n的矩阵。A的伴随矩阵可按如下步骤定义：

定义：A关于第i行第j列的余子式(记作M_ij)是去掉A的第i行第j列之后得到的(n-1)*(n-1)矩阵的行列式。

定义：A关于第i行第j列的代数余子式是：C_ij=(-1)^i+jM_ij。

定义：A的余子矩阵是一个n*n的矩阵C，使得其第i行第j列的元素是A关于第i行第j列的代数余子式。

引入以上的概念后，可以定义：矩阵A的伴随矩阵是A的余子矩阵的转置矩阵：adj(A)= C^T.

也就是说，A的伴随矩阵是一个n*n的矩阵(记作adj(A))，使得其第i行第j列的元素是A关于第j行第i列的代数余子式：[adj(A)]_ij=C_ji。

伴随矩阵性质：对n*n的矩阵A和B，有：

(1)、adj(I) = I;

(2)、adj(AB) = adj(B)adj(A);

(3)、adj(A^T) = adj(A)^T;

(4)、det(adj(A)) = det(A)^n-1;

(5)、adj(kA) = k^n-1adj(A);

(6)、当n>2时，adj(adj(A)) = (detA)^n-2A;

(7)、如果A可逆，那么adj(A^-1) = adj(A)^-1 = A / detA;

(8)、如果A是对称矩阵，那么其伴随矩阵也是对称矩阵；如果A是反对称矩阵，那么当n为偶数时，A的伴随矩阵也是反对称矩阵，n为奇数时则是对称矩阵；

(9)、如果A是(半)正定矩阵，那么其伴随矩阵也是(半)正定矩阵；

(10)、如果矩阵A和B相似，那么adj(A)和adj(B)也相似；

(11)、如果n>2，那么非零矩阵A是正交矩阵当且仅当adj(A) = ±A^T.

以上内容整理自 维基百科

以下是用C++实现的求伴随矩阵：

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

// 计算行列式
template<typename _Tp>
_Tp determinant(const std::vector<std::vector<_Tp>>& mat, int N)
{
	if (mat.size() != N) {
		fprintf(stderr, "mat must be square matrix\n");
		return -1;
	}
	for (int i = 0; i < mat.size(); ++i) {
		if (mat[i].size() != N) {
			fprintf(stderr, "mat must be square matrix\n");
			return -1;
		}
	}

	_Tp ret{ 0 };

	if (N == 1) return mat[0][0];

	if (N == 2) {
		return (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
	}
	else {
		// first col
		for (int i = 0; i < N; ++i) {
			std::vector<std::vector<_Tp>> m(N - 1);
			std::vector<int> m_rows;
			for (int t = 0; t < N; ++t) {
				if (i != t) m_rows.push_back(t);
			}
			for (int x = 0; x < N - 1; ++x) {
				m[x].resize(N - 1);
				for (int y = 0; y < N - 1; ++y) {
					m[x][y] = mat[m_rows[x]][y + 1];
				}
			}
			int sign = (int)pow(-1, 1 + i + 1);
			ret += mat[i][0] * sign * determinant<_Tp>(m, N - 1);
		}
	}

	return ret;
}

// 计算伴随矩阵
template<typename _Tp>
int adjoint(const std::vector<std::vector<_Tp>>& mat, std::vector<std::vector<_Tp>>& adj, int N)
{
	if (mat.size() != N) {
		fprintf(stderr, "mat must be square matrix\n");
		return -1;
	}
	for (int i = 0; i < mat.size(); ++i) {
		if (mat[i].size() != N) {
			fprintf(stderr, "mat must be square matrix\n");
			return -1;
		}
	}

	adj.resize(N);
	for (int i = 0; i < N; ++i) {
		adj[i].resize(N);
	}

	for (int y = 0; y < N; ++y) {
		std::vector<int> m_cols;
		for (int i = 0; i < N; ++i) {
			if (i != y) m_cols.push_back(i);
		}

		for (int x = 0; x < N; ++x) {
			std::vector<int> m_rows;
			for (int i = 0; i < N; ++i) {
				if (i != x) m_rows.push_back(i);
			}

			std::vector<std::vector<_Tp>> m(N - 1);
			for (int i = 0; i < N - 1; ++i) {
				m[i].resize(N - 1);
			}
			for (int j = 0; j < N - 1; ++j) {
				for (int i = 0; i < N - 1; ++i) {
					m[j][i] = mat[m_rows[j]][m_cols[i]];
				}
			}

			int sign = (int)pow(-1, x + y);
			adj[y][x] = sign * determinant<_Tp>(m, N-1);
		}
	}

	return 0;
}

template<typename _Tp>
void print_matrix(const std::vector<std::vector<_Tp>>& mat)
{
	int rows = mat.size();
	for (int y = 0; y < rows; ++y) {
		for (int x = 0; x < mat[y].size(); ++x) {
			fprintf(stderr, "  %f  ", mat[y][x]);
		}
		fprintf(stderr, "\n");
	}
	fprintf(stderr, "\n");
}

int test_adjoint_matrix()
{
	std::vector<float> vec{5, -2, 2, 7, 1, 0, 0, 3, -3, 1, 5, 0, 3, -1, -9, 4 };
	const int N{ 4 };
	if (vec.size() != (int)pow(N, 2)) {
		fprintf(stderr, "vec must be N^2\n");
		return -1;
	}

	std::vector<std::vector<float>> arr(N);
	for (int i = 0; i < N; ++i) {
		arr[i].resize(N);

		for (int j = 0; j < N; ++j) {
			arr[i][j] = vec[i * N + j];
		}
	}

	std::vector<std::vector<float>> adj;
	int ret = adjoint<float>(arr, adj, N);

	fprintf(stderr, "source matrix: \n");
	print_matrix<float>(arr);
	fprintf(stderr, "adjoint matrx: \n");
	print_matrix<float>(adj);

	return 0;
}

输出结果如下：

GitHub：https://github.com/fengbingchun/NN_Test