当前位置:   article > 正文

六种语言求解特征值和特征向量(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 [ a11a12...a1na21a22...a2n............an1an2...ann \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 (λEA)x =0是包含n个未知数n个方程的齐次线性方程组,相当于存在不全为0的n个数使得矩阵 ( λ E − A ) (\lambda E -A) (λEA)的列向量线性相关,那么矩阵 ( λ E − A ) (\lambda E -A) (λEA)非满秩,即 r a n k ( λ E − A ) < n rank(\lambda E -A)<n rank(λEA)<n        另一方面,线性代数中齐次线性方程组存在非零解充要条件是:系数矩阵行列式的值为零。因此得到下面的系数行列式:
∣ λ E − A ∣ = 0 |\lambda E -A|=0 λEA=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 [ λa11a12...a1na21λa22...a2n............an1an2...λann \right ]=\lambda^{n}+\alpha_{1}\lambda^{n-1}+\alpha_{2}\lambda^{n-2}+...+\alpha_{n-1}\lambda+\alpha_{n}=0 λEA=λa11a21...an1a12λa22...an2............a1na2n...λann=λn+α1λn1+α2λn2+...+αn1λ+α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=1naii λ 1 λ 2 ⋯ λ n = ∣ A ∣ \lambda_{1}\lambda_{2}\cdots\lambda_{n}=|A| λ1λ2λn=A

2、六种编程语言求解特征值和特征向量代码

        为了求解如下方阵A特征值特征向量,下面分别利用MatlabPythonC++JavaC#、和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 [ 3451037179262485247194127832653123169574362128495812181617345229279931579083 \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);
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26

        Matlab代码运行结果如下图所示:
在这里插入图片描述

2.2 Python代码及运行结果(Numpy)

        这里利用Python 3.6numpy依赖库来对矩阵进行特征值和特征向量运算,具体代码如下。

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=' ')
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28

        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;
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47

        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();
        }
	}
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59

        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();
            }
        }
    }
}
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53

        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)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9

        R代码运行结果如下图所示:
在这里插入图片描述

3、总结

        总的来说,MatlabPythonR同为脚本语言,可解释性强,易于执行,易于验证,直观高效;而JavaC#C++都为面向对象的高级编程语言,需要依托主函数才能执行代码,效率较高,计算速度快,但往往需要额外依托第三方库或者自己开发代码来执行矩阵运算。
        同时实际测试发现:六种编程语言对特征值的计算都十分准确,但特征向量结果不同,具体不同之处在于:MatlabPython计算的特征向量结果相同,为复数JavaC#计算的特征向量结果相同,是实数C++计算的特征向量结果为复数,但与MatlabPython的计算结果略有差异;R处理结果相差较大。最终利用Matlab软件验证特征值特征向量的结果如下图所示。因此,无论利用何种编程语言进行计算,如果对计算结果的精度要求较高,建议与专业软件进行结果比对和相互验证,从而保障计算的准确性,这也告诉我们第三方库也无法保证算法的高度准确,因为它们所采用的计算方式可能有所差异

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/笔触狂放9/article/detail/78155
推荐阅读
相关标签
  

闽ICP备14008679号