特征值和特征向量求解(2022.4.29)
引言
在使用层次分析法(AHP)来确定指标权重时,需要计算二维矩阵的特征值和特征向量,同时特征值和特征向量在计算机图像处理
、物理学
等领域都要十分重要的应用,下面首先简单介绍其定义
和性质
,之后利用六种编程语言进行实例求解。
1、矩阵的特征值和特征向量
A
=
[
a
11
a
12
.
.
.
a
1
n
a
21
a
22
.
.
.
a
2
n
.
.
.
.
.
.
.
.
.
.
.
.
a
n
1
a
n
2
.
.
.
a
n
n
]
A=\left [ \begin{matrix} &a_{11} &a_{12} &... &a_{1n}\\ &a_{21} &a_{22} &... &a_{2n}\\ &... &... &... &...\\ &a_{n1} &a_{n2} &... &a_{nn}\\ \end{matrix} \right ]
A=⎣⎢⎢⎡a11a21...an1a12a22...an2............a1na2n...ann⎦⎥⎥⎤ 对于n
阶方阵A
,若存在数λ
和一个n
维的非零列向量
x
⃗
\vec{x}
x,满足
A
x
⃗
=
λ
x
⃗
A\vec{x}=\lambda \vec{x}
Ax=λx,则称
x
⃗
\vec{x}
x为矩阵A
的对应λ
特征值的特征向量,E
为单位矩阵(主对角线元素均为1,其余为0)。
一方面,
(
λ
E
−
A
)
x
⃗
=
0
(\lambda E -A) \vec{x}=0
(λE−A)x=0是包含n
个未知数n
个方程的齐次线性方程组,相当于存在不全为0的n
个数使得矩阵
(
λ
E
−
A
)
(\lambda E -A)
(λE−A)的列向量线性相关,那么矩阵
(
λ
E
−
A
)
(\lambda E -A)
(λE−A)非满秩,即
r
a
n
k
(
λ
E
−
A
)
<
n
rank(\lambda E -A)<n
rank(λE−A)<n 另一方面,线性代数中齐次线性方程组存在非零解的充要条件是:系数矩阵行列式的值为零。因此得到下面的系数行列式:
∣
λ
E
−
A
∣
=
0
|\lambda E -A|=0
∣λE−A∣=0
∣
λ
E
−
A
∣
=
[
λ
−
a
11
−
a
12
.
.
.
−
a
1
n
−
a
21
λ
−
a
22
.
.
.
−
a
2
n
.
.
.
.
.
.
.
.
.
.
.
.
−
a
n
1
−
a
n
2
.
.
.
λ
−
a
n
n
]
=
λ
n
+
α
1
λ
n
−
1
+
α
2
λ
n
−
2
+
.
.
.
+
α
n
−
1
λ
+
α
n
=
0
|\lambda E -A|=\left [ \begin{matrix} &\lambda -a_{11} &-a_{12} &... &-a_{1n}\\ &-a_{21} &\lambda -a_{22} &... &-a_{2n}\\ &... &... &... &...\\ &-a_{n1} &-a_{n2} &... &\lambda -a_{nn}\\ \end{matrix} \right ]=\lambda^{n}+\alpha_{1}\lambda^{n-1}+\alpha_{2}\lambda^{n-2}+...+\alpha_{n-1}\lambda+\alpha_{n}=0
∣λE−A∣=⎣⎢⎢⎡λ−a11−a21...−an1−a12λ−a22...−an2............−a1n−a2n...λ−ann⎦⎥⎥⎤=λn+α1λn−1+α2λn−2+...+αn−1λ+αn=0 上式也称为方阵A
的特征多项式,此多项式是一个关于参数λ
的n
次多项式。在复数域内(实数和虚数),n
次代数方程有且仅有n
个解。那么特征多项式的n个解可设为
λ
1
、
λ
2
、
⋯
、
λ
n
\lambda_{1}、\lambda_{2}、\cdots、\lambda_{n}
λ1、λ2、⋯、λn,这些解也称特征根,它们具有如下性质:
(
λ
−
λ
1
)
(
λ
−
λ
2
)
⋯
(
λ
−
λ
n
)
=
0
(\lambda -\lambda_{1})(\lambda -\lambda_{2})\cdots (\lambda -\lambda_{n})=0
(λ−λ1)(λ−λ2)⋯(λ−λn)=0
λ
1
+
λ
2
+
.
.
.
+
λ
n
=
∑
i
=
1
n
a
i
i
\lambda_{1}+\lambda_{2}+...+\lambda_{n}=\sum_{i=1}^{n}a_{ii}
λ1+λ2+...+λn=i=1∑naii
λ
1
λ
2
⋯
λ
n
=
∣
A
∣
\lambda_{1}\lambda_{2}\cdots\lambda_{n}=|A|
λ1λ2⋯λn=∣A∣
2、六种编程语言求解特征值和特征向量代码
为了求解如下方阵A
的特征值和特征向量,下面分别利用Matlab
、Python
、C++
、Java
、C#
、和R
这六种编程语言进行代码实现,同时对比测试结果。
A
=
[
3
4
5
10
3
7
1
7
9
2
6
2
4
8
5
2
4
7
1
9
4
12
78
32
65
31
23
16
95
74
36
21
28
49
58
12
18
16
17
34
52
29
2
7
99
31
57
90
83
]
A=\left [ \begin{matrix} &3 &4 &5 &10 &3 &7 &1\\ &7 &9 &2 &6 &2 &4 &8\\ &5 &2 &4 &7 &1 &9 &4\\ &12 &78 &32 &65 &31 &23 &16\\ &95 &74 &36 &21 &28 &49 &58\\ &12 &18 &16 &17 &34 &52 &29\\ &2 &7 &99 &31 &57 &90 &83\\ \end{matrix} \right ]
A=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡3751295122492787418752432361699106765211731321312834577492349529018416582983⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
2.1 Matlab代码及运行结果
作为工业界和学术界的神器,MathWorks
旗下的Matlab备受科研人员、工程师和开发者的喜爱,功能强大且有着“矩阵实验室”称号,在数学运算
、图像处理
、信号仿真
方面有着其他软件无法比拟的计算优势,尽管有时其效率不是太高,但其拥有着毋庸置疑的准确率和精度。
n = 7;
A=[3,4,5,10,3,7,1;
7,9,2,6,2,4,8;
5,2,4,7,1,9,4;
12,78,32,65,31,23,16;
95,74,36,21,28,49,58;
12,18,16,17,34,52,29;
2,7,99,31,57,90,83];
[eigenvector,x] = eig(A);
lambda = diag(x);
fprintf('特征值:');
for i=1:n
if i == n
fprintf(['λ',num2str(i),'= ',num2str(lambda(i)),'\n']);
else
fprintf([' ','λ',num2str(i),'= ',num2str(lambda(i)),' ']);
end
end
for i=1:n
if i == n
fprintf(['特征向量',num2str(i),'\n']);
else
fprintf([' ','特征向量',num2str(i),' ']);
end
end
disp(eigenvector);
Matlab
代码运行结果如下图所示:
2.2 Python代码及运行结果(Numpy)
这里利用Python 3.6
和numpy
依赖库来对矩阵进行特征值和特征向量运算,具体代码如下。
import numpy as np
n = 7
A = np.array([[3,4,5,10,3,7,1],
[7,9,2,6,2,4,8],
[5,2,4,7,1,9,4],
[12,78,32,65,31,23,16],
[95,74,36,21,28,49,58],
[12,18,16,17,34,52,29],
[2,7,99,31,57,90,83]])
print('方阵A:\n{}'.format(A))
eigenvalue, eigenvector = np.linalg.eig(A)
print('特征值', end='')
for i in range(0, n):
if i == n-1:
print('λ{}={:.4f}'.format(i + 1, eigenvalue[i]))
else:
print('λ{}={:.4f}'.format(i+1, eigenvalue[i]), end=', ')
for i in range(0, n):
if i == n-1:
print(' 特征向量{} '.format(i + 1))
else:
print(' 特征向量{} '.format(i + 1), end=' ')
for i in range(0, n):
for j in range(0, n):
if j == n-1:
print('{:.4f} '.format(eigenvector[i][j]))
else:
print('{:.4f} '.format(eigenvector[i][j]), end=' ')
Python
代码运行结果如下图所示:
2.3 C++代码及运行结果(Eigen)
Eigen作为线性代数相关计算的C++
模板库广受开发人员的喜爱,下面使用Eigen
库来计算矩阵特征值和特征向量,下载eigen-3.4.0
后,在VS 2015
中将其目录下的Eigen文件夹添加到项目属性的库目录中即可使用。
#include <iostream>
#include <iomanip>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
const int n = 7;
Matrix<double, n, n> A = MatrixXd::Random(n, n);
double a[] = { 3,4,5,10,3,7,1,
7,9,2,6,2,4,8,
5,2,4,7,1,9,4,
12,78,32,65,31,23,16,
95,74,36,21,28,49,58,
12,18,16,17,34,52,29,
2,7,99,31,57,90,83};
A = Map<Matrix<double, n, n, RowMajor> >(a);
std::cout << "方阵A:" << std::endl;
std::cout << A << std::endl;
EigenSolver<Matrix<double, n, n>> eigensolver(A);
MatrixXcd evecs = eigensolver.eigenvectors();//矩阵特征向量complex复数矩阵
MatrixXcd evals = eigensolver.eigenvalues();
//cout << evals << endl;
//cout << evecs << endl;
std::cout << "特征值:";
for (int i = 0; i < n; i++)
if(evals(i).imag() >= 0)
cout << "λ" << i + 1 << "="<< setprecision(4) <<evals(i).real() <<"+"<< setprecision(4)<< evals(i).imag()<< "i ";
else
cout << "λ" << i + 1 << "="<< setprecision(4) << evals(i).real() << setprecision(4) << evals(i).imag() << "i ";
std::cout << std::endl;
for(int i=0;i<n;i++)
std::cout << "特征向量" << (i+1) <<" ";
std::cout << std::endl;
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
if(evecs(i, j).imag() >= 0)
std::cout << setprecision(4) << evecs(i, j).real() <<"+"<< setprecision(4) << evecs(i, j).imag()<<"i ";
else
std::cout << setprecision(4) << evecs(i, j).real() << setprecision(4) << evecs(i, j).imag() << "i ";
}
std::cout << std::endl;
}
return 0;
}
C++
代码运行结果如下图所示:
2.4 Java代码及运行结果(ujmp)
这里使用第三方jar
包:ujmp-core-0.3.0.jar
来实现Java代码,jar包下载地址为ujmp-core-0.3.0.jar,之后在Java
项目中需要将其添加到Build Path
中,才能在代码中导入。
package com.test;
import java.text.DecimalFormat;
import org.ujmp.core.DenseMatrix;
import org.ujmp.core.Matrix;
import org.ujmp.core.doublematrix.impl.ArrayDenseDoubleMatrix2D;
public class EigenValue_EigenVector
{
public static void main(String[] args)
{
int n = 7;
double[][] a = {{3,4,5,10,3,7,1},
{7,9,2,6,2,4,8},
{5,2,4,7,1,9,4},
{12,78,32,65,31,23,16},
{95,74,36,21,28,49,58},
{12,18,16,17,34,52,29},
{2,7,99,31,57,90,83}};
Matrix A = DenseMatrix.Factory.importFromArray(a);
System.out.print("A:");
System.out.println(A);
Matrix[] eigenValueDecomposition = A.eig();
Matrix eigen_vector = eigenValueDecomposition[0];
Matrix eigen_value = eigenValueDecomposition[1];
double[][] eig_val= eigen_value.toDoubleArray();
double[][] eig_vec= eigen_vector.toDoubleArray();
DecimalFormat decimalFormat = new DecimalFormat("###################.####");
System.out.print("特征值");
for(int i=0;i<n;i++)
{
String val = String.valueOf(decimalFormat.format(eig_val[i][i]));
for(int j=0;j<n;j++)
{
if(j !=i && eig_val[j][i] != 0)
if(eig_val[j][i] >= 0)
val += "+"+String.valueOf(decimalFormat.format(eig_val[j][i]))+"i";
else
val += String.valueOf(decimalFormat.format(eig_val[j][i]))+"i";
}
System.out.print("λ"+Integer.toString(i+1)+"= "+val+" ");
}
System.out.println();
for(int i=0;i<n;i++)
{
System.out.print(" 特征向量"+Integer.toString(i+1)+" ");
}
System.out.println();
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
{
System.out.print("\t"+String.valueOf(decimalFormat.format(eig_vec[i][j]))+"\t");
}
System.out.println();
}
}
}
Java
代码运行结果如下图所示:
2.5 C#代码及运行结果(MathNet)
这里利用了适用于C#
语言的经典数学dll
库:MathNet.Numerics.dll
,在编写代码时需要将其添加到引用当中,以便在代码中调用其接口。
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Numerics;
using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;//必须
using MathNet.Numerics.LinearAlgebra.Complex;
using MathNet.Numerics.LinearAlgebra.Factorization;
namespace SolveMatrixEigen
{
class Program
{
public static void Main(string[] args)
{
Console.WriteLine("Hello,World!");
const int n = 7;
double[,] a = new double[n, n]{{3,4,5,10,3,7,1},
{7,9,2,6,2,4,8},
{5,2,4,7,1,9,4},
{12,78,32,65,31,23,16},
{95,74,36,21,28,49,58},
{12,18,16,17,34,52,29},
{2,7,99,31,57,90,83}};
MathNet.Numerics.LinearAlgebra.Matrix<double> A = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix.OfArray(a);
MathNet.Numerics.LinearAlgebra.Factorization.Evd<double> eigen = A.Evd();
Vector<Complex> Eigen_Value = eigen.EigenValues;//矩阵的特征值
MathNet.Numerics.LinearAlgebra.Matrix<double> Eigen_Vector = eigen.EigenVectors;//矩阵的特征向量
Console.Write("特征值:");
for (int i = 0; i < n; i++)
{
string eigenvalue_i = "";
if(Eigen_Value[i].Imaginary >= 0)
eigenvalue_i = String.Format("{0:F4}+{1:F4}i", Eigen_Value[i].Real, Eigen_Value[i].Imaginary);
else
eigenvalue_i = String.Format("{0:F4}{1:F4}i", Eigen_Value[i].Real, Eigen_Value[i].Imaginary);
Console.Write("λ" + (i + 1).ToString() + "=" + eigenvalue_i+" ");
}
Console.WriteLine();
for (int i = 0; i < n; i++)
Console.Write(" 特征向量{0} ", i + 1);
Console.WriteLine();
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
Console.Write(String.Format(" {0:F4} ",Eigen_Vector[i,j]));
Console.WriteLine();
}
}
}
}
C#
代码运行结果如下图所示:
2.6 R代码及运行结果
Windows下的R安装包下载,安装成功后出现RGui
启动程序,进入后可直接编写脚本代码。
a <- c(3,4,5,10,3,7,1,7,9,2,6,2,4,8,5,2,4,7,1,9,4,12,78,32,65,31,23,16,95,74,36,21,28,49,58,12,18,16,17,34,52,29,2,7,99,31,57,90,83)
A <- matrix(a,7,7)
print("方阵A")
print(A)
ev <- eigen(A)
print("特征值:")
print(ev$values)
print("特征向量:")
print(ev$vectors)
R
代码运行结果如下图所示:
3、总结
总的来说,Matlab
、Python
和R
同为脚本语言,可解释性强,易于执行,易于验证,直观高效;而Java
、C#
和C++
都为面向对象的高级编程语言,需要依托主函数才能执行代码,效率较高,计算速度快,但往往需要额外依托第三方库或者自己开发代码来执行矩阵运算。
同时实际测试发现:六种编程语言对特征值的计算都十分准确,但特征向量结果不同,具体不同之处在于:Matlab
和Python
计算的特征向量结果相同,为复数;Java
和C#
计算的特征向量结果相同,是实数;C++
计算的特征向量结果为复数,但与Matlab
和Python
的计算结果略有差异;R
处理结果相差较大。最终利用Matlab
软件验证特征值和特征向量的结果如下图所示。因此,无论利用何种编程语言进行计算,如果对计算结果的精度要求较高,建议与专业软件进行结果比对和相互验证,从而保障计算的准确性,这也告诉我们第三方库也无法保证算法的高度准确,因为它们所采用的计算方式可能有所差异。
![]() | ![]() |
![]() | ![]() |
![]() | ![]() |
![]() |