【数值计算方法】矩阵特征值与特征向量的计算(一):Jacobi 旋转法及其Python实现

文章目录

  • 一、Jacobi 旋转法
    • 1. 基本思想
    • 2. 计算过程演示
    • 3. 注意事项
  • 二、Python实现
    • 迭代过程(调试)

  矩阵的特征值(eigenvalue)和特征向量(eigenvector)在很多应用中都具有重要的数学和物理意义。Jacobi 旋转法是一种用于计算对称矩阵特征值和特征向量的迭代方法。

  本文将详细介绍 Jacobi 旋转法的基本原理和步骤,通过一个具体的矩阵示例演示其应用过程,并给出其Python实现。

一、Jacobi 旋转法

  Jacobi 旋转法的每一次迭代中,需要选择一个非对角元素最大的位置,然后构造相应的旋转矩阵,进行相似变换,使得矩阵逐渐对角化。

  • 对称矩阵是一个实数矩阵,其转置与自身相等。
  • 对于一个方阵 A A A,如果存在标量 λ λ λ 和非零向量 v v v,使得 A v = λ v Av = λv Av=λv,那么 λ λ λ 就是 A A A 的特征值, v v v 就是对应于 λ λ λ 的特征向量。

1. 基本思想

  Jacobi 旋转法的基本思想是通过一系列的相似变换,逐步将对称矩阵对角化,使得非对角元素趋于零。这个过程中,特征值逐渐浮现在对角线上,而相应的特征向量也被逐步找到。下面是 Jacobi 旋转法的基本步骤:

  1. 选择旋转角度: 选择一个旋转角度 θ,通常使得旋转矩阵中的非对角元素为零,从而实现对角化,通常选择非对角元素中绝对值最大的那个作为旋转的目标。

  2. 构造旋转矩阵: 构造一个旋转矩阵 J,该矩阵为单位矩阵,只有对应于选择的非对角元素的位置上有两个非零元素,其余位置上为零。这两个非零元素的值由旋转角度 θ 决定,例如,对于 2x2 矩阵,旋转矩阵可以表示为:
    J = [ cos ⁡ ( θ ) − sin ⁡ ( θ ) sin ⁡ ( θ ) cos ⁡ ( θ ) ] J = \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix} J=[cos(θ)sin(θ)sin(θ)cos(θ)]

  3. 相似变换: 计算相似变换矩阵 P P P,即 P T A P P^TAP PTAP,其中 A A A 是原始矩阵, P P P 是旋转矩阵,计算过程如下:

P T A P = [ cos ⁡ ( θ ) sin ⁡ ( θ ) − sin ⁡ ( θ ) cos ⁡ ( θ ) ] T [ a 11 a 12 a 12 a 22 ] [ cos ⁡ ( θ ) − sin ⁡ ( θ ) sin ⁡ ( θ ) cos ⁡ ( θ ) ] P^TAP = \begin{bmatrix} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{bmatrix}^T \begin{bmatrix} a_{11} & a_{12} \\ a_{12} & a_{22} \end{bmatrix} \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix} PTAP=[cos(θ)sin(θ)sin(θ)cos(θ)]T[a11a12a12a22][cos(θ)sin(θ)sin(θ)cos(θ)]

  通过矩阵相乘计算,我们可以得到 P T A P P^TAP PTAP 中的非对角元素,假设这两个元素分别位于矩阵的 (1,2) 和 (2,1) 的位置。令 a i j a_{ij} aij 为这两个元素,即 a i j = a 12 = a 21 a_{ij}= a_{12} = a_{21} aij=a12=a21

  接下来,我们希望通过选择合适的 θ \theta θ使得 a i j a_{ij} aij 变为零,从而达到对角化的目的,即 a 12 = a 21 a_{12} = a_{21} a12=a21,进一步可推导出

θ = 1 2 arctan ⁡ ( 2 ⋅ a i j a i i − a j j ) \theta = \frac{1}{2} \arctan\left(\frac{2 \cdot a_{ij}}{a_{ii} - a_{jj}}\right) θ=21arctan(aiiajj2aij)

  • a i i = a j j a_{ii}=a_{jj} aii=ajj,则使用 a r c c o t arccot arccot形式
  1. 迭代: 重复步骤 1-3,直到矩阵 A 的非对角元素都趋于零或满足一定的精度要求。

  2. 提取特征值和特征向量: 对角线上的元素即为矩阵 A 的特征值,而 P 中的列向量即为对应于这些特征值的特征向量。

2. 计算过程演示

  对于矩阵
A = [ 2 − 1 0 − 1 2 − 1 0 − 1 2 ] A = \begin{bmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{bmatrix} A= 210121012

  我们首先找到非对角元素中绝对值最大的元素,这里我们以 (2,1) 为例,计算旋转角度和旋转矩阵。

  1. 选择旋转角度:

      计算旋转角度 θ \theta θ公式:
    θ = 1 2 arctan ⁡ ( 2 ⋅ a i j a i i − a j j ) \theta = \frac{1}{2} \arctan\left(\frac{2 \cdot a_{ij}}{a_{ii} - a_{jj}}\right) θ=21arctan(aiiajj2aij)其中, a i i a_{ii} aii a j j a_{jj} ajj 分别是矩阵的对角元素,而 a i j a_{ij} aij 是非对角元素,即 a 21 a_{21} a21。 在这个例子中, a 21 = − 1 a_{21} = -1 a21=1 a 11 = a 22 = 2 a_{11} = a_{22} = 2 a11=a22=2

    θ = 1 2 arctan ⁡ ( − 2 0 ) = − π 4 \theta = \frac{1}{2} \arctan\left(\frac{-2}{0}\right) = -\frac{\pi}{4} θ=21arctan(02)=4π

  2. 构造旋转矩阵:

    构造旋转矩阵 ( J ):

J = [ cos ⁡ ( θ ) − sin ⁡ ( θ ) sin ⁡ ( θ ) cos ⁡ ( θ ) ] J = \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix} J=[cos(θ)sin(θ)sin(θ)cos(θ)]

对于 θ = − π 4 \theta = -\frac{\pi}{4} θ=4π

J = [ cos ⁡ ( − π 4 ) − sin ⁡ ( − π 4 ) sin ⁡ ( − π 4 ) cos ⁡ ( − π 4 ) ] J = \begin{bmatrix} \cos\left(-\frac{\pi}{4}\right) & -\sin\left(-\frac{\pi}{4}\right) \\ \sin\left(-\frac{\pi}{4}\right) & \cos\left(-\frac{\pi}{4}\right) \end{bmatrix} J=[cos(4π)sin(4π)sin(4π)cos(4π)]

计算得:

J = [ 2 2 2 2 − 2 2 2 2 ] J = \begin{bmatrix} \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} \\ -\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} \end{bmatrix} J=[22 22 22 22 ]

  1. 相似变换:

    计算相似变换矩阵 P P P

    P T A P P^T A P PTAP

    在这里, P P P就是构造的旋转矩阵 J J J

  2. 迭代:

    重复上述步骤,直到矩阵足够接近对角矩阵。

  这个过程会一步步地使矩阵趋近于对角矩阵,对角线上的元素就是矩阵的特征值,而相应的列向量就是对应的特征向量。由于计算较为繁琐,我在这里只展示了一个迭代的过程。在实际应用中,你需要进行多次迭代,直到满足精度的要求。
在这里插入图片描述
在这里插入图片描述

3. 注意事项

  Jacobi 旋转法的优点是可以用于任意大小的对称矩阵,但其缺点是迭代次数较多,计算量较大。在实际应用中,通常会结合其他方法来提高计算效率。

二、Python实现

import numpy as npdef jacobi_rotation(A):n = A.shape[0]tolerance = 1e-10max_iterations = 1000eigenvectors = np.eye(n)for _ in range(max_iterations):# 寻找最大的非对角元素max_off_diag = np.max(np.abs(np.triu(A, k=1)))if max_off_diag < tolerance:break  # 达到收敛条件# 找到最大元素的索引indices = np.unravel_index(np.argmax(np.abs(np.triu(A, k=1))), A.shape)i, j = indices# 计算旋转角度theta = 0.5 * np.arctan2(2 * A[i, j], A[i, i] - A[j, j])# 构造旋转矩阵J = np.eye(n)J[i, i] = J[j, j] = np.cos(theta)J[i, j] = -np.sin(theta)J[j, i] = np.sin(theta)# 执行相似变换A = np.dot(np.dot(J.T, A), J)# 更新特征向量eigenvectors = np.dot(eigenvectors, J)# 提取特征值eigenvalues = np.diag(A)return eigenvalues, eigenvectors# 示例矩阵
A = np.array([[2, -1, 0],[-1, 2, -1],[0, -1, 2]])# 执行 Jacobi 旋转
eigenvalues, eigenvectors = jacobi_rotation(A)print("特征值:", eigenvalues)
print("特征向量:")
np.set_printoptions(precision=4, suppress=True)
print(eigenvectors)

在这里插入图片描述

迭代过程(调试)

  • 第一次:
    在这里插入图片描述
  • 第二次:在这里插入图片描述
    ………
  • 第九次:
    在这里插入图片描述

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/155573.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

软件工程理论与实践 (吕云翔) 第十二章 软件测试概述课后习题及其答案解析

第十二章 软件测试概述 1.判断题 &#xff08;1&#xff09;软件测试是对软件规格说明、软件设计和编码的最全面也是最后的审查。&#xff08;) 解析&#xff1a;软件测试是一个独立的活动&#xff0c;旨在评估软件的质量和发现潜在的问题&#xff0c;而不仅仅是对规格说明、…

六、流量监管、流量整形

流量监管、流量整形 流量监管、流量整形1.1.定义1.2.简单流分类应用场景举例1.3.简单流分类的配置1.4.复杂流分类1.5.复杂流分类在产品中的实现 ————————————————— 流量监管、流量整形 当报文的发送速率大于接收速率&#xff0c;或者下游设备的接口速率小于…

使用Python合并pdf文件

import PyPDF2 # 创建一个空的 PDF 文档对象 def pdf_merge(path_a,path_b,path_save):merged_pdf PyPDF2.PdfMerger()# 打开第一个 PDF 文件&#xff0c;读取其内容并追加到新创建的 PDF 文档对象中with open(path_a, rb) as file1:merged_pdf.append(file1)# 打开第二个 PDF…

GitHub 2023报告-开源和AI的现状

GitHub 2023报告-开源和AI的现状 深入探讨人工智能如何与开源互动&#xff0c;以及未来几年可能出现的趋势。 背景介绍 2023年&#xff0c;开源已成为全球软件开发的标准。无论是大公司还是小团队&#xff0c;都广泛使用开源技术进行项目开发。此外&#xff0c;随着机器学习和…

Docker部署FLASK Unicorn并配置Nginx

1. 安装相关依赖 flask3.0.0 pymysql1.1.0 #我自己需要的 flask_cors4.0.0 gunicorn21.2.0 gevent23.9.12. 配置Gunicorn 新建gunicorn.conf.py bind 0.0.0.0:5418 # 绑定的IP地址和端口 workers 8 # 同时执行的进程数&#xff0c;推荐为当前CPU个数*21 worker_class&qu…

【Java】java | CacheManager | redisCacheManager

一、说明 1、查询增加缓存&#xff0c;使用Cacheable注解 2、项目中已经用到了ehcache&#xff0c;现在需求是两个都用 二、备份配置 1、redisConfig增加代码 Bean("redisCacheManage")Primarypublic CacheManager redisCacheManager(RedisConnectionFactory fact…

项目集成的关键过程总结

项目集成是确保项目各部分协调一致&#xff0c;共同实现项目目标的过程。以下是项目集成的关键过程&#xff1a; 制定集成计划 制定集成计划是项目集成的第一步。该计划应明确集成的目标、范围、时间表、资源分配和风险管理策略。在制定集成计划时&#xff0c;需要考虑项目干…

CKD TransBTS:用于脑肿瘤分割的具有模态相关交叉注意的临床知识驱动混合转换器

CKD-TransBTS: Clinical Knowledge-Driven Hybrid Transformer With Modality-Correlated Cross-Attention for Brain Tumor Segmentation CKD TransBTS&#xff1a;用于脑肿瘤分割的具有模态相关交叉注意的临床知识驱动混合转换器背景贡献实验方法how radiologists diagnose b…

【19年扬大真题】已知a数组int a[ ]={1,2,3,4,5,6,7,8,9,10},编写程序,求a数组中偶数的个数和偶数的平均值

【18年扬大真题】 已知a数组int a[ ]{1,2,3,4,5,6,7,8,9,10}&#xff0c;编写程序&#xff0c;求a数组中偶数的个数和偶数的平均值 int main() {int arr[10] { 1,2,3,4,5,6,7,8,9,10 };int os 0;//偶数个数int sum 0;//偶数和float ave 0;//偶数平均值for (int i 0;i <…

关于2023年11月25日PMI认证考试准考信下载及考场规定等事项通知

各位考生&#xff1a;为保证参加2023年11月25日PMI项目管理资格认证考试的每位考生都能顺利进入考场参加考试&#xff0c;请完整阅读本通知内容。 一、关于准考信下载为确保您顺利进入考场参加11月份考试&#xff0c;请及时登录本网站个人系统下载并打印准考信&#xff0c;准考…

Google codelab WebGPU入门教程源码<6> - 使用计算着色器实现计算元胞自动机之生命游戏模拟过程(源码)

对应的教程文章: https://codelabs.developers.google.com/your-first-webgpu-app?hlzh-cn#7 对应的源码执行效果: 对应的教程源码: 此处源码和教程本身提供的部分代码可能存在一点差异。点击画面&#xff0c;切换效果。 class Color4 {r: number;g: number;b: number;a…

论文阅读:“基于特征检测与深度特征描述的点云粗对齐算法”

文章目录 摘要简介相关工作粗对齐传统的粗对齐算法基于深度学习的粗对齐算法 特征检测及描述符构建 本文算法ISS 特征检测RANSAC 算法3DMatch 算法 实验结果参考文献 摘要 点云对齐是点云数据处理的重要步骤之一&#xff0c;粗对齐则是其中的难点。近年来&#xff0c;基于深度…

【机器学习Python实战】线性回归

&#x1f680;个人主页&#xff1a;为梦而生~ 关注我一起学习吧&#xff01; &#x1f4a1;专栏&#xff1a;机器学习python实战 欢迎订阅&#xff01;后面的内容会越来越有意思~ ⭐内容说明&#xff1a;本专栏主要针对机器学习专栏的基础内容进行python的实现&#xff0c;部分…

C++ Qt 学习(十):Qt 其他技巧

1. 带参数启动外部进程 QProcess 用于启动外部进程int QProcess::execute(const QString &program, const QStringList &arguments);QObject *parent; ... QString program "./path/to/Qt/examples/widgets/analogclock"; QStringList arguments; argument…

当前系统并无桌面环境,或无显示器,无法显示远程桌面,您需要自行安装X11桌面环境,或者使用终端文件功能

ToDesk远程遇到的问题如上图&#xff0c;换向日葵直接黑屏&#xff1b; 问题原因 截止发文时间&#xff0c;Todesk只支持X11协议&#xff0c;没有适配最新的Wayland协议&#xff0c;所以我们需要把窗口系统调整为X11才可以。 解决方法 修改配置文件&#xff0c;关闭wayland su…

Vellum|SOP —— Vellum Constraints

目录 约束类型 Distance Along Edges&#xff08;distance&#xff09; Bend Across Triangles&#xff08;bend&#xff09; Cloth&#xff08;distance、bend&#xff09; Hair&#xff08;stretchshear、bendtwist&#xff09; String&#xff08;distance、angle&…

生成指定位数强Lucas校验伪素数-Arnault1995构造法

Arnault在1995年的论文《Constructing Carmichael Numbers which are Strong Pseudoprimes to Several Bases》提出了一种用于构造强Lucas校验伪素数的方法&#xff0c;本文将对其方法做具体的实现分析。 文章目录 1.Lucas素数测试1.1 Lucas序列1.2 Lucas定理1.3 Lucas素数测试…

Spring Cloud Netflix微服务组件-Hystrix

目录 Hystrix的主要功能 传统容错手段 超时机制 应用容错三板斧 超时机制 舱壁隔离 熔断降级 侵入式Command用法 改进版一&#xff1a;ribbon与hystrix组合 改进版二&#xff1a;feign与hystrix组合 Hystrix三态转换图 源码分析 流程图 原理流程图 核心实现流程…

MySQL 的执行原理(一)

5.1 单表访问之索引合并 我们前边说过 MySQL 在一般情况下执行一个查询时最多只会用到单个二级 索引&#xff0c;但存在有特殊情况&#xff0c;在这些特殊情况下也可能在一个查询中使用到多个二 级索引&#xff0c;MySQL 中这种使用到多个索引来完成一次查询的执行方法称之为&…

基于Android个人理财家庭收支系统uniAPP+vue -hbuiderx-微信小程序vj9n8

摘 要 个人理财APP设计的目的是为用户提供对活动信息和提醒信息管理等方面的平台。 与PC端应用程序相比&#xff0c;个人理财的设计旨在为用户提供一个个人理财信息的管理。用户可以通过APP及时收入信息和支出信息&#xff0c;并对提醒消息进行查看等。 个人理财APP是在Androi…