【计算方法与科学建模】矩阵特征值与特征向量的计算(一):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 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)

在这里插入图片描述

迭代过程(调试)

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

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

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

缓存雪崩、击穿、穿透_解决方案

文章目录 缓存雪崩、击穿、穿透1.缓存雪崩造成缓存雪崩解决缓存雪崩 2. 缓存击穿造成缓存击穿解决缓存击穿 3.缓存穿透造成缓存穿透解决缓存穿透 缓存雪崩、击穿、穿透 一般用户数据存储于磁盘&#xff0c;读写速度慢。 使用redis作为缓存&#xff0c;相当于数据缓存在内存&a…

ILI9225 TFT显示屏16位并口方式驱动

所用屏及资料如后图&#xff1a; ILI9225&#xff0c;176*220&#xff0c;8位或16位并口屏&#xff0c;IM0接GND&#xff0c;电源及背光接3.3v 主控&#xff1a;CH32V307驱动&#xff08;库文件和STM32基本一样&#xff09; 一、源码 ILI9225.c #include "ILI9225.h&quo…

二分查找——经典题目合集

文章目录 &#x1f99c;69. x 的平方根&#x1f33c;题目&#x1f33b;算法原理&#x1f337;代码实现 &#x1f433;35. 搜索插入位置&#x1f33c;题目&#x1f33b;算法原理&#x1f337;代码实现 &#x1f9ad;852. 山脉数组的峰顶索引&#x1f33c;题目&#x1f33b;算法原…

IBM ELM—系统工程全生命周期管理平台

产品概述 Engineering Lifecycle Management是IBM提供的工程全生命周期管理组合工具&#xff0c;帮助企业降低开发成本&#xff0c;应对开发挑战并更快地发展其流程和实践。 随着产品变得更加复杂且数字化&#xff0c;传统的工程开发不再能及时且有效地满足系统工程的复杂度&a…

AD7021C 触摸感应加灯光调节芯片IC 可用于触摸台灯、触摸玩具灯等

AD7021C触摸感应 IC 是为实现人体触摸界面而设计的集成电路。可替代机械式轻触按键&#xff0c;实现防水防尘、密封隔离、坚固美观的操作界面。使用该芯片可以实现 LED 灯光亮度调节&#xff0c;方案所需的外围电路简单&#xff0c;操作方便。确定好灵敏度选择电容&#xff…

最新PHP熊猫头图片表情斗图生成源码

这是一款能生成熊猫头表情斗图的自适应系统源码&#xff0c;无论是在电脑还是手机上都可以正常使用&#xff01;这个源码集成了搜狗搜索图片接口&#xff0c;可以轻松地一键搜索数百万张图片&#xff0c;并且还包含了表情制作等功能模块。对于一些新站来说&#xff0c;这是一个…

商务俄语学习,柯桥基础入门教学,千万别小看俄语中的“что”

1、что до (чего) 至于 例&#xff1a; что до меня, то я не могу согласиться 至于我&#xff0c;我不能同意。 А что до зимовки... Ты приедешь в этом году? 说到冬天和过冬…你今年回来吗…

【unity实战】unity3D中的PRG库存系统和换装系统(附项目源码)

文章目录 先来看看最终效果前言素材简单绘制库存UI前往mixamo获取人物模型动画获取一些自己喜欢的装备物品模型库存系统换装系统装备偏移问题添加消耗品最终效果源码完结 先来看看最终效果 前言 之前2d的换装和库存系统我们都做过不少了&#xff0c;这次就来学习一个3d版本的&…

利用Python进行数据分析【送书第六期:文末送书】

&#x1f468;‍&#x1f393;博主简介 &#x1f3c5;云计算领域优质创作者   &#x1f3c5;华为云开发者社区专家博主   &#x1f3c5;阿里云开发者社区专家博主 &#x1f48a;交流社区&#xff1a;运维交流社区 欢迎大家的加入&#xff01; &#x1f40b; 希望大家多多支…

MQTT协议消息代理服务远程连接

目录 1. Linux 搭建 Mosquitto 2. Linux 安装Cpolar 3. 创建MQTT服务公网连接地址 4. 客户端远程连接MQTT服务 5. 代码调用MQTT服务 6. 固定连接TCP公网地址 7. 固定地址连接测试 Mosquitto是一个开源的消息代理&#xff0c;它实现了MQTT协议版本3.1和3.1.1。它可以在不…

HCIP-三、VRRP+BFD

三、VRRPBFD 实验拓扑实验需求及解法 实验拓扑 实验需求及解法 //本实验模拟某公司网关冗余结构&#xff0c;按以下要求完成配置&#xff1a; //1.如图所示&#xff0c;配置 R1/2/3 的设备名称及 IP 地址。 //2.内外网通信。 //2.1 在 R1/2 上配置默认路由&#xff0c;保证 R1…

爆款文章有诀窍,内容创作者如何能持续产出优质内容

内容营销人有没有这么一种共鸣&#xff1a;10 万 那么多&#xff0c;为什么不能多我一个&#xff1f; 通常&#xff0c;我们把浏览量 / 阅读量高、转评赞数量高的内容看作爆款&#xff0c;而数据如果达到 10 万 则是超级爆款。因为&#xff0c;阅读量高意味着内容得到了大量的曝…

2023年03月 Scratch(二级)真题解析#中国电子学会#全国青少年软件编程等级考试

Scratch等级考试(1~4级)全部真题・点这里 一、单选题(共25题,每题2分,共50分) 第1题 小猫的程序如图所示,积木块的颜色与球的颜色一致。点击绿旗执行程序后,下列说法正确的是?( ) A:小猫一直在左右移动,嘴里一直说着“抓到了”。 B:小猫会碰到球,然后停止。…

[点云分割] 条件欧氏聚类分割

介绍 条件欧氏聚类分割是一种基于欧氏距离和条件限制的点云分割方法。它通过计算点云中点与点之间的欧氏距离&#xff0c;并结合一定的条件限制来将点云分割成不同的区域或聚类。 在条件欧氏聚类分割中&#xff0c;通常会定义以下两个条件来判断两个点是否属于同一个聚类&…

记录小白第一次EDUsrc:任意密码漏洞

目录 一、漏洞说明&#xff1a; 二、漏洞复现&#xff1a; 三、漏洞修复建议&#xff1a; 一、漏洞说明&#xff1a; xxxx学院身份认证系统有严重的逻辑设计缺陷&#xff1a;账户登录、手机登录、密码找回三个接口找到n个逻辑漏洞包括任意账号密码修改、信息泄露&#xff0…

AIGC ChatGPT4总结SQL优化细节操作

数据库SQL优化是一个复杂的过程,它通常涉及到许多不同的技术和方法。以下是一些常用的SQL优化策略: 1. **索引使用**:索引可以极大地加速查询速度。但是,索引并不总是有好处的,因为它们需要额外的空间来存储,并且在插入和更新数据时可能会减慢速度。因此,选择正确的字段…

windows系统安装ubuntu22.04虚拟机

镜像文件准备 镜像文件 官网 企业开源和Linux | Ubuntu 镜像下载地址 https://cn.ubuntu.com/download/server/step1 选择合适的版本下载 虚拟机安装 文件-- 新建虚拟机 选择镜像 修改安装路径 修改大小&#xff0c;最好60g&#xff0c;大一点 设置用户信息 设置虚拟机网络…

羊大师提示,羊奶都有哪些惊人功效?

羊奶不仅是一种美味的健康饮品&#xff0c;在近年来备受瞩目的的健康圈子里&#xff0c;羊奶还被赋予了更多的功效&#xff0c;成为一种备受推崇的保健品。羊奶不但富含营养&#xff0c;而且还有着非常多的益处&#xff0c;它能够用来美容、保健&#xff0c;甚至还可以治疗某些…

Java8新特性 ----- Lambda表达式和方法引用/构造器引用详解

前言 在讲一下内容之前,我们需要引入函数式接口的概念 什么是函数式接口呢? 函数式接口&#xff1a;有且仅有一个抽象方法的接口 java中函数式编程的体现就是Lambda表达式,你可以认为函数式接口就是适用于Lambda表达式的接口. 也可以加上注解来在编译层次上限制函数式接口 Fun…

音频采集的相关基础知识

本文引注: https://zhuanlan.zhihu.com/p/652629744 1.麦克风的种类 (1)模拟麦克风 ECM麦克风&#xff1a;驻极体电容麦克风(ECM)&#xff0c;典型的汽车ECM麦克风是一种将ECM单元与小型放大器电路整合在单个外壳中的装置。放大器提供一个模拟信号&#xff0c;其电压电平允许…