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

文章目录

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

在这里插入图片描述

迭代过程(调试)

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

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

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

相关文章

06、基于内容的过滤算法Tensorflow实现

06、基于内容的过滤算法Tensorflow实现 开始学习机器学习啦&#xff0c;已经把吴恩达的课全部刷完了&#xff0c;现在开始熟悉一下复现代码。对这个手写数字实部比较感兴趣&#xff0c;作为入门的素材非常合适。 05、基于梯度下降的协同过滤算法中已经介绍了协同过滤算法的基…

用最少数量的箭引爆气球[中等]

优质博文&#xff1a;IT-BLOG-CN 一、题目 有一些球形气球贴在一堵用XY平面表示的墙面上。墙面上的气球记录在整数数组points&#xff0c;其中points[i] [xstart, xend]表示水平直径在xstart和xend之间的气球。你不知道气球的确切y坐标。一支弓箭可以沿着x轴从不同点完全垂直…

Springboot-注册注解【springboot常用注解】

1.组件注册 1.1 使用的注解 Configuration:普通配置类,替代以前的配置文件,配置类本身也是容器的组件|SpringBootConfiguration:Springboot配置类,与Configuration功能一样|Bean:替代以前的Bean标签,如果没有在Bean标签内定义名字,则默认组件的名字为方法名,可以直接修改注解…

在gitlab上使用server_hooks

文章目录 1. 前置条件2. Git Hook2.1 Git Hook 分为两部分&#xff1a;本地和远程2.1.1 本地 Git Hook&#xff0c;由提交和合并等操作触发&#xff1a;2.1.2 远程 Git Hook&#xff0c;运行在网络操作上&#xff0c;例如接收推送的提交&#xff1a; 3. 操作步骤3.1 对所有的仓…

Linux下的文件IO之系统IO

1. 知识点 读入写出&#xff0c;切记以我们程序为中心向文件或者别的什么东西读入写出&#xff08;输入流输出流&#xff09; 人话就是 文件向我们程序就是读入 程序向文件或者别的什么就是写出 2. open打开文件 open.c /****************************************************…

C++ 学习之函数成员指针的一个小细节

看看下面的代码&#xff0c;你能看出错误吗 class A { public:void fun(){}}; int main() {A a;void (A:: * p)() &A::fun;(*p)(); } 这段代码在调用成员函数时存在问题。正确的方式是使用对象来调用成员函数&#xff0c;而不是通过指针。以下是修正后的代码&#xff1a…

STM32CubeIDE(CUBE-MX hal库)----定时器

系列文章目录 STM32CubeIDE(CUBE-MX hal库)----初尝点亮小灯 STM32CubeIDE(CUBE-MX hal库)----按键控制 STM32CubeIDE(CUBE-MX hal库)----串口通信 文章目录 系列文章目录前言一、定时器二、使用步骤三、HAL库实验代码三、标准库代码 前言 STM32定时器是一种多功能外设&#…

异常 Exception 02

异常 Exception 02 六、异常处理1、基本介绍2、异常处理的方式3、示意图 try-catchthrows1、介绍2、注意事项 自定义异常1、基本概念2、自定义异常的步骤3、实例4、throw和throws的区别 六、异常处理 1、基本介绍 异常处理就是当异常发生时,对异常处理的方式。 2、异常处理的…

以STM32CubeMX创建DSP库工程方法一

以STM32CubeMX创建DSP库工程方法 略过时钟树的分配和UART的创建等&#xff0c;直接进入主题生成工程文件 它们中的文件功能如下&#xff1a; 1&#xff09;BasicMathFunctions 基本数学函数&#xff1a;提供浮点数的各种基本运算函数&#xff0c;如向量加减乘除等运算。 2&…

VBA代码解决方案第8讲:用FindPrevious进行重复搜索及利用LIKE查找

《VBA代码解决方案》(版权10028096)这套教程是我最早推出的教程&#xff0c;目前已经是第三版修订了。这套教程定位于入门后的提高&#xff0c;在学习这套教程过程中&#xff0c;侧重点是要理解及掌握我的“积木编程”思想。要灵活运用教程中的实例像搭积木一样把自己喜欢的代码…

货代FOB条款卖方必备的知识:发货人都要承担哪些费用呢?

据统计&#xff0c;中国出口中以FOB成交的占到70%&#xff0c;但专家指出&#xff1a;FOB对出口商的风险更大&#xff0c;有可能造成货、款两空的结局。 目前我国出口合同以FOB价格条款成交的比例越来越大&#xff0c;而且收货人指定船公司的少&#xff0c;指定境外货代的多&am…

3款厉害的小工具,小黑子都在用!

大家好&#xff0c;我是 Javapub。 程序员与普通人最大的区别是什么&#xff0c;当然是会使用工具。基于一些同学经常问我的问题&#xff0c;接下来给大家分享几款我经常使用的工具&#xff0c;主打一个提升效率。 第一款 Everything 用 windwos 的同学都体会过&#xff0c;…

带大家做一个,易上手的家常土豆片

还是先从冰箱里那一块猪瘦肉 搞一点蒜和生姜 切成小块 装进一个碗里 这里一点就够了 一条绿皮辣椒 切片 三个左右干辣椒 随便切两刀 让它小一点就好了 一起装一个碗 一大一小两个土豆切片 猪肉切片 起锅烧油 然后 下肉翻炒 等肉变颜色捞出来 然后放入土豆 和小半碗水 让…

JeecgBoot低代码开发—Vue3版前端入门教程

JeecgBoot低代码开发—Vue3版前端入门教程 后端接口配置VUE3 必备知识1.vue3新特性a. https://v3.cn.vuejs.org/b.setup的用法c.ref 和 reactive 的用法d.新版 v-model 的用法e.script setup的用法 2.TypeScript基础 后端接口配置 如何修改后台项目路径 http://127.168.3.52:8…

【玩转 EdgeOne】| 腾讯云下一代边缘加速CDN EdgeOne 是安全加速界的未来吗?

目录 前言边缘加速与安全加固边缘计算与CDN的融合EdgeOne优秀的安全特性EdgeOne卓越的性能表现灵活的配置和管理生态系统的支持与发展技术创新与未来展望EdgeOne试用结束语 前言 在当下互联网的迅猛发展的时刻&#xff0c;云计算和边缘计算技术的快速发展为网络加速领域带来了…

awk从放弃到入门(10):awk内置函数

awk从放弃到入门&#xff08;10&#xff09;&#xff1a;awk内置函数 算数函数字符串函数其他函数 本博文转载自 这篇文章中的知识点是建立在前文的基础上的&#xff0c;如果你还没有掌握前文中的知识&#xff0c;请先参考之前的文章。 注&#xff1a;在阅读这篇文章之前&…

Abbyy FineReader16最新版本有哪些新功能?

在数字化时代&#xff0c;数据处理和转换变得非常重要&#xff0c;Abbyy FineReader 就是一款专门用于处理、转换和识别图像和 PDF 文件的软件。在本文中&#xff0c;我们将会详细介绍 Abbyy FineReader 的功能以及适合使用该软件的电脑。 ABBYY Finereader 16-安装包下载如下&…

滴滴2023.11.27P0级故障技术复盘回顾(k8s的的错?)

本文从滴滴官方恢复及技术公众号带大家从技术角度复盘这次事故 目录 1. 背景 2. 滴滴官方消息 3. 问题分析及定位 4.网传的k8s及解析 5.k8s引发的思考&#xff1a;举一反三&#xff0c;怎么避免再次出现 6.近段时间其他平台崩溃回顾 1. 背景 11 月 27 晚约 10 点&#xf…

TFIDF、BM25、编辑距离、倒排索引

TFIDF TF刻画了词语t对某篇文档的重要性&#xff0c;IDF刻画了词语t对整个文档集的重要性

最简单的链路追踪收集器

链路追踪可帮助您快速了解程序服务之间的调用关系&#xff0c;并快速洞悉内部发生的情况。主流的链路追踪系统有zipkin,jaeger,skywalking等&#xff0c;由于opentelemetry的存在&#xff0c;都具有opentelemetry的转换器。 我们利用opentelemetry来进行zipkin,jaeger,skywalk…