【数值计算方法(黄明游)】矩阵特征值与特征向量的计算(三):Jacobi 旋转法【理论到程序】

news2024/11/15 13:57:00

文章目录

  • 一、Jacobi 旋转法
    • 1. 基本思想
    • 2. 计算过程演示
  • 二、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 = [ 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. 迭代:

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

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

在这里插入图片描述
在这里插入图片描述

二、Python实现

import numpy as np


def jacobi_rotation(A):
    n = A.shape[0]
    tolerance = 1e-10
    max_iterations = 1000
    eigenvectors = 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.coloradmin.cn/o/1273103.html

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

相关文章

华为云之云桌面Workspace的使用体验

华为云之云桌面Workspace的使用体验 一、云桌面Workspace介绍1.云桌面简介2.云桌面特点3. 云桌面应用场景①远程移动办公②协同办公③安全办公④公用终端⑤图形制作渲染 二、本次实践介绍1. 本次实践目的2. 本次实践环境 三、购买云桌面1. 进入华为云的云桌面购买界面2. 选择购…

智能优化算法应用:基于生物地理学算法无线传感器网络(WSN)覆盖优化 - 附代码

智能优化算法应用&#xff1a;基于生物地理学算法无线传感器网络(WSN)覆盖优化 - 附代码 文章目录 智能优化算法应用&#xff1a;基于生物地理学算法无线传感器网络(WSN)覆盖优化 - 附代码1.无线传感网络节点模型2.覆盖数学模型及分析3.生物地理学算法4.实验参数设定5.算法结果…

【蓝桥杯选拔赛真题70】Scratch输入输出数字 少儿编程scratch图形化编程 蓝桥杯创意编程选拔赛真题解析

目录 scratch输入输出数字 一、题目要求 编程实现 二、案例分析 1、角色分析

文件夹重命名技巧:用关键词替换文件夹名称指定内容的右侧文字

在日常生活中&#xff0c;经常要管理大量的文件夹&#xff0c;这时候掌握一些文件夹重命名的技巧就非常实用。例如文件夹重命名时&#xff0c;经常要将一些通用的文字替换成其他关键词&#xff0c;以便更好地标识和分类文件夹。而用关键词替换文件夹名称指定内容的右侧文字&…

智能优化算法应用:基于纵横交叉算法无线传感器网络(WSN)覆盖优化 - 附代码

智能优化算法应用&#xff1a;基于纵横交叉算法无线传感器网络(WSN)覆盖优化 - 附代码 文章目录 智能优化算法应用&#xff1a;基于纵横交叉算法无线传感器网络(WSN)覆盖优化 - 附代码1.无线传感网络节点模型2.覆盖数学模型及分析3.纵横交叉算法4.实验参数设定5.算法结果6.参考…

亚马逊云科技re:Invent Peter DeSantis演讲,数据规模拓展无极限引领Serverless构建之路

re:lnvent 2023 Peter DeSantis主题演讲&#xff0c;数据规模拓展无极限引领Serverless构建之路&#xff08;Road to Serverless&#xff09;。 Logical Qubit全新发布&#xff1a;量子计算硬件&#xff0c;6倍的量子纠错效率提升。 Amazon全新发布Redshift Serverless&#xf…

RedisTemplate中使用scan方法代替keys指令

keys * 这个命令千万别在生产环境乱用&#xff0c;危险级别不亚于flushdb。因为Redis是单线程操作&#xff0c;在数据特别庞大的情况下。Keys会引发Redis锁&#xff08;数据过多一直查询处理&#xff09;&#xff0c;并且增加Redis的CPU占用。很多公司的运维都是禁止了这个命令…

Python实现FA萤火虫优化算法优化循环神经网络回归模型(LSTM回归算法)项目实战

说明&#xff1a;这是一个机器学习实战项目&#xff08;附带数据代码文档视频讲解&#xff09;&#xff0c;如需数据代码文档视频讲解可以直接到文章最后获取。 1.项目背景 萤火虫算法&#xff08;Fire-fly algorithm&#xff0c;FA&#xff09;由剑桥大学Yang于2009年提出 , …

狗都会配的SNAT和DNAT配置

1 SNAT 1.1 SNAT SNAT原理与应用:. SNAT 应用环境:局域网主机共享单个公网IP地址接入Internet (私有IP不能在Internet中正常路由) SNAT原理:源地址转换&#xff0c;根据指定条件修改数据包的源IP地址&#xff0c;通常被叫做源映谢 SNAT转换前提条件: 1.局域网各主机已正确设…

13-Vue基础之自定义指令与插槽的使用

个人名片&#xff1a; &#x1f60a;作者简介&#xff1a;一名大二在校生 &#x1f921; 个人主页&#xff1a;坠入暮云间x &#x1f43c;座右铭&#xff1a;懒惰受到的惩罚不仅仅是自己的失败&#xff0c;还有别人的成功。 &#x1f385;**学习目标: 坚持每一次的学习打卡 文章…

Apipost推出IDEA插件,代码写完直接调试

IDEA是一款功能强大的集成开发环境&#xff08;IDE&#xff09;&#xff0c;它可以帮助开发人员更加高效地编写、调试和部署软件应用程序。我们在编写完接口代码后需要进行接口调试等操作&#xff0c;一般需要打开额外的调试工具。 今天给大家介绍一款IDEA插件&#xff1a;Api…

mysql 日志分析

程序启动标志 可以直接全局搜索&#xff0c;查看启动了几次 可以看到总共11次&#xff0c;当前是第2次 如何判断mysql是正常关闭&#xff0c;手动启动的 下图中启动之前出现 Shutdown complete打印说明启动之前是正常关闭的

2024 本田 CBR1000RR-R Fireblade SP

本田在米兰车展期间也发布了自家的旗舰仿赛&#xff0c;虽然不是重大改款&#xff0c;只是更新一些内容&#xff0c;性能方面有所小提升&#xff0c;但是毕竟人家是本田&#xff0c;毕竟人家是火刃&#xff0c;还是要尊重一下&#xff0c;写一篇内容的。 新款还是有一些亮点&am…

京东秒杀之秒杀实现

1 登录判断 用户在未登录状态下可以查看商品列别以及秒杀商品详情&#xff0c;但不可以在未登录状态进行秒杀商品的操作&#xff0c;当用户点击开始秒杀时&#xff0c;进行登陆验证 <!DOCTYPE html> <head><title>商品详情</title><meta http-eq…

C++ -- 每日选择题 -- Day2

第一题 1. 下面代码中sizeof(A)结果为&#xff08;&#xff09; #pragma pack(2) class A {int i;union U{char str[13];int i;}u;void func() {};typedef char* cp;enum{red,green,blue}color; }; A&#xff1a;20 B&#xff1a;21 C&#xff1a;22 D&#xff1a;24 答案及解析…

[NOIP2016 普及组] 回文日期

枚举好题&#xff0c;直接枚举答案 看看在不在范围内就行了 注意二月份 92200229是合法的~ 82200228也是合法的&#xff01; #include<bits/stdc.h> using namespace std;map<int,int>mp;int main() {mp[1] mp[3] mp[5] mp[7] mp[8] mp[10] mp[12] 31;mp[…

kafka C++实现消费者

文章目录 1 Kafka 消费者的逻辑2 Kafka 的C API2.1 RdKafka::Conf2.2 RdKafka::Event2.3 RdKafka::EventCb2.4 RdKafka::TopicPartition2.5 RdKafka::RebalanceCb2.6 RdKafka::Message2.7 RdKafka::KafkaConsumer&#xff08;核心&#xff09; 3 Kafka 消费者客户端开发3.1 必要…

为大家收集了一些最常用的Python包

我们从最常用的 Python 包入手&#xff0c;去解答上述这个问题。最初&#xff0c;我列出过去一年在 PyPI 上下载次数最多的 Python 包。接下来&#xff0c;深入研究其用途、它们之间的关系和它们备受欢迎的原因。 1、Urllib3 下载次数&#xff1a;8.93 亿 Urllib3是一个 Pyt…

MongoDB快速入门及其SpringBoot实战

MongoDB快速入门及其SpringBoot实战 MongoDB简介 MongoDB 是一个基于分布式文件存储的数据库。由 C 语言编写。旨在为 WEB 应用提供可扩展的高性能数据存储解决方案。 MongoDB是一个开源、高性能、无模式的文档型数据库&#xff0c;当初的设计就是用于简化开发和方便扩展&am…

numpy模块安装方法

https://www.bilibili.com/video/BV1qN411R7V2/?spm_id_from333.337.search-card.all.click&vd_sourcefb8dcae0aee3f1aab700c21099045395