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

news2024/10/5 23:30:29

文章目录

  • 一、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)

在这里插入图片描述

迭代过程(调试)

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

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

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

相关文章

Linux 环境配置小白入门

Linux从 全栈开发centOS 7 到 运维 一 Linux 入门概述1.1 操作系统1.2 Linux 简介1.3 Linux 系统组成1.4 Linux 发行版1.5 Linux 应用领域1.6 Linux vs Windows 二 虚拟机2.1 虚拟机介绍2.2 VMware WorkStation 安装2.3 VMware WorkStation 配置检查2.3 安装 CentOS 72.3.1 安装…

Git 远程仓库(Github)

目录 添加远程库 查看当前的远程库 提取远程仓库 推送到远程仓库 删除远程仓库 Git 并不像 SVN 那样有个中心服务器。 目前我们使用到的 Git 命令都是在本地执行&#xff0c;如果你想通过 Git 分享你的代码或者与其他开发人员合作。 你就需要将数据放到一台其他开发人员…

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

Scratch等级考试(1~4级)全部真题・点这里 一、单选题(共25题,每题2分,共50分) 第1题 小猫在沙漠中旅行好不容易找到了一杯水,初始位置如下图所示,下面哪个程序可以帮助它成功喝到水? A: B: C: D:

LR学习笔记——基本面板

文章目录 面板介绍色彩调整区域明暗调整区域纹理及质感色彩饱和 面板介绍 面板如上图所示 基本可分为几个板块&#xff1a;色彩、明暗、纹理及质感、色彩饱和 色彩调整区域 色温&#xff1a;由蓝色和黄色控制色调&#xff1a;由绿色和洋红控制 互补色&#xff1a;蓝色对黄色&…

opencv-形态学处理

通过阈值化分割可以得到二值图&#xff0c;但往往会出现图像中物体形态不完整&#xff0c;变的残缺&#xff0c;可以通过形态学处理&#xff0c;使其变得丰满&#xff0c;或者去除掉多余的像素。常用的形态学处理算法包括&#xff1a;腐蚀&#xff0c;膨胀&#xff0c;开运算&a…

Altium Designer学习笔记7

PCB封装库的制作&#xff1a; 距离的测量&#xff1a; 各个焊盘的位置&#xff1a; 直插元件选择Multi-Layer。如果贴片元件的则选择顶层Top-Layer&#xff0c;或者Bottom-Layer。 形状是方形&#xff0c;尺寸是2mm*2mm。 孔的尺寸是1.4mm。 则该器件就制作完成。 TSSOP28封装…

Java精品项目源码基于SpringBoot的樱花短视频平台(v66)

Java精品项目源码基于SpringBoot的樱花短视频平台(v66) 大家好&#xff0c;小辰今天给大家介绍一个樱花短视频平台&#xff0c;演示视频公众号&#xff08;小辰哥的Java&#xff09;对号查询观看即可 文章目录 Java精品项目源码基于SpringBoot的樱花短视频平台(v66)难度指数&…

Python3,必备数据可视化之:数据交互可视化

数据可视化之交互可视化 1、引言2、交互可视化介绍2.1 Bokeh2.1.1 基本定义2.1.2 常用功能2.1.3 安装2.1.4 代码示例 2.2 Plotly2.2.1 基本定义2.1.2 常用功能2.1.3 安装2.2.4 代码示例 2.3 Bokeh与Plotly 差异点 3、总结 1、引言 小屌丝&#xff1a;鱼哥&#xff0c;我发现一…

详解Python安装requests库的实例代码

文章目录 前言基本用法基本的get请求带参数的GET请求解析json关于Python技术储备一、Python所有方向的学习路线二、Python基础学习视频三、精品Python学习书籍四、Python工具包项目源码合集①Python工具包②Python实战案例③Python小游戏源码五、面试资料六、Python兼职渠道 前…

腾讯云COS+picgo+typora 图床搭建与自动上传

1、腾讯云 COS 腾讯云活动 COS新用户专享 COS 操作步骤 1、点击 创建桶&#xff0c;完善信息 点击下一步&#xff0c;剩下的配置可自己配置 2、picgo 官网地址 2.3.1版本下载地址 现在稳定版本是2.3.1 相关连接 腾讯云密钥设置地址picgo官网地址2.3.1版本下载地址

4-11 四个数排序

#include<stdio.h> int main(){int t,a,b,c,d;printf("请输入四个数&#xff1a;");scanf("%d %d %d %d",&a,&b,&c,&d);printf("a%d,b%d,c%d,d%d\n",a,b,c,d);if(a>b){ta;ab;bt;}if(a>c){ta;ac;ct;}if(a>d){ta;a…

注册亚马逊美国买家号需要些什么资料?

注册亚马逊美国买家号需要准备邮箱、美国手机号、美国地址及能支付的支付卡。准备好之后进入亚马逊美国站进行点击注册&#xff0c;按照格式填写好之后即可注册成功了。 而如果想要注册大量买家号&#xff0c;可以使用亚马逊鲲鹏系统进行自动化操作&#xff0c;想要自动化更顺畅…

【STM32】TF卡FTA32文件系统

一、SD卡介绍 1.SD简介 本质&#xff1a;NandFlash控制芯片 2.SD卡存储容量等级 3.FAT文件系统的使用 4.SD卡速度等级 5.SD卡驱动方式 1.SDIO&&SD 1&#xff09;SDIO接口通信线&#xff1a;CLK/CMD/DAT0-3&#xff08;数据传输线4根&#xff09; 2&#xff09;SPI接口…

这样做出来的电子杂志好看又精美,不信你也来试试!

大家会不会有一种困惑&#xff0c;为什么别人制作的电子杂志那么精美高级&#xff0c;能翻页&#xff0c;能分享到微信、微博等&#xff0c;而自己制作的电子杂志却是平平无奇&#xff1f; 在这个快节奏的时代&#xff0c;人们的生活方式也发生了翻天覆地的变化。而使用FLBOOK…

(swjtu西南交大)数据库实验(数据库需求分析):音乐软件数据管理系统

实验内容&#xff1a; 数据库需求分析&#xff1a;各用户组需求描述&#xff0c;绘出数据流图&#xff08;详细案例参见教材p333~p337&#xff0c;陶宏才&#xff0c;数据库原理及设计&#xff0c;第三版&#xff09;&#xff1b; 一、选题背景 近年来&#xff0c;“听歌”逐…

Python-函数传参与数据类型

Python中&#xff0c;函数参数传递是通过对象的引用进行的&#xff0c;我们可以进行下面的验证。 def use_name(val):print("name id :%s" % (id(val)))val "hanshu1"print("name id modified :%s" % (id(val)))def test_ref():name "ha…

CyNix

CyNix 一、主机发现和端口扫描 主机发现&#xff0c;靶机地址192.168.80.146 arp-scan -l端口扫描&#xff0c;只开放了80和6688端口 nmap -A -p- -sV 192.168.80.146二、信息收集 访问80端口 路径扫描 gobuster dir -u http://192.168.80.146/ -w /usr/share/wordlists/dir…

吴恩达《机器学习》9-4-9-6:实现注意:展开参数、梯度检验、随机初始化

一、实现注意:展开参数 在上一个视频中&#xff0c;讨论了使用反向传播算法计算代价函数的导数。在本视频中&#xff0c;将简要介绍一个实现细节&#xff0c;即如何将参数从矩阵展开为向量。这样做是为了在高级最优化步骤中更方便地使用这些参数。 二、梯度检验 在神经网络中…

如何制作动态表情包?一个方法快学起来

在当代的通讯工具中&#xff0c;动态表情包已经是人们日常交流不可缺少的一部分了。但是&#xff0c;很多时候网络上常见的动态表情包不能够很好表达出我们的需求时应该怎么办呢&#xff1f;这时候&#xff0c;我们可以使用gif动图制作&#xff08;https://www.gif.cn/&#xf…

深入理解强化学习——马尔可夫决策过程:马尔可夫决策过程和马尔可夫过程/马尔可夫奖励过程的区别

分类目录&#xff1a;《深入理解强化学习》总目录 《深入理解强化学习——马尔可夫决策过程》系列前面的文章讨论到的马尔可夫过程和马尔可夫奖励过程都是自发改变的随机过程&#xff0c;而如果有一个外界的“刺激”来共同改变这个随机过程&#xff0c;就有了马尔可夫决策过程&…