【数值计算方法(黄明游)】矩阵特征值与特征向量的计算(三):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 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/186725.shtml

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

相关文章

LeetCode Hot100 438.找到字符串中所有字母异位词

题目&#xff1a; 给定两个字符串 s 和 p&#xff0c;找到 s 中所有 p 的 异位词 的子串&#xff0c;返回这些子串的起始索引。不考虑答案输出的顺序。 异位词 指由相同字母重排列形成的字符串&#xff08;包括相同的字符串&#xff09;。 代码&#xff1a; class Solution …

位图和布隆过滤器(C++)

位图和布隆过滤器 一、位图1. 引入2. 概念3. 代码实现setreset完整代码 4. 位图的应用 二、布隆过滤器1. 引入2. 概念3. 逻辑结构4. 特点5. 代码实现6. 布隆过滤器的应用 三、哈希切割 一、位图 1. 引入 当面对海量数据需要处理时&#xff0c;内存不足以加载这些数据&#xf…

The Sandbox 携手 Sandsoft,与 Nuqtah 合作推动沙特阿拉伯的 Web3 发展

新的合作伙伴关系将增强创作者的能力&#xff0c;促进区块链生态系统的包容性。 The Sandbox 及其合作伙伴 Sandsoft 是移动游戏开发商和发行商&#xff0c;也是 AAA 人才驱动的投资者&#xff0c;他们非常高兴地宣布与 Nuqtah 建立新的合作伙伴关系&#xff0c;Nuqtah 是中东和…

numpy模块安装方法

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

Linux:Ubuntu系统安装软件

本次以安装vim为例 sudo apt-get remove vim //卸载vim sudo apt-get install vim //安装vim sudo apt-cache show vim //获取vim软件信息安装时间较长。 安装完成后&#xff0c;执行下第三条指令&#xff0c;测试下是否安装成功即可。

在gazebo里搭建一个livox mid360 + 惯导仿真平台测试 FAST-LIO2

在gazebo里搭建一个livox mid360 惯导仿真平台测试 FAST-LIO2 前言立方体平台加入 livox mid360 激光雷达加入IMU模块调整底盘大小 并设计调用接口测试 Fast-Lio2 前言 livox mid360 在官网一直没有货&#xff0c;在gazebo里可以仿真该雷达形式的点云。 但是其只发布雷达的数…

Spire.Office 8.11.2 for NET fix Crack

内容摘自来自互联网------或者SDK官方本身手册 Spire.Doc for .NET A professional Word .NET library designed to create, read, write, convert and print Word document files in any .NET ( C#, VB.NET, ASP.NET, .NET Core, Xamarin ) application with fast and high qu…

Aurora8B10B(一) 从IP配置界面学习Aurora

一. 简介 哈喽&#xff0c;大家好&#xff0c;好久没有给大家写FPGA技术的文章&#xff0c;是不是已经忘记我是做FPGA的啦&#xff0c;O(∩_∩)O哈哈~。 这里将会给大家分享我学习到的第一个高速接口Aurora8B10B&#xff0c;有点复杂&#xff0c;但不是特别复杂&#xff0c;对…

使用vscode的remotessh插件远程连接的时候被要求重复输入密码

问题描述&#xff1a; 需要远程连接服务器&#xff0c;使用ssh&#xff0c;我用到的是vscode里面的remotessh插件。配置好config以后 HostHostNameUserPortIdentifyFile进入到了vscode的密码登录界面&#xff0c;但是一直被要求循环输入密码&#xff0c;很奇怪&#xff0c;去…

论文阅读——DINOv

首先是关于给了提示然后做分割的一些方法的总结&#xff1a; 左边一列是prompt类型&#xff0c;右边一列是使用各个类型的prompt的模型。这些模型有分为两大类&#xff1a;Generic和Refer&#xff0c;通用分割和参考分割。Generic seg 是分割和提示语义概念一样的所有的物体&am…

LLM之Agent(二):BabyAGI的详细教程

BabyAGI是一个 AI 支持的任务管理系统&#xff08;Python脚本&#xff09;&#xff0c;使用 OpenAI 和 Pinecone API 创建, 优先级排序和执行任务。该系统背后的主要思想是基于先前任务的结果和预定义的目标创建任务。脚本然后使用 OpenAI 的自然语言处理&#xff08;NLP&#…

leetCode 93.复原 IP 地址 + 回溯算法 + 图解 + 笔记

93. 复原 IP 地址 - 力扣&#xff08;LeetCode&#xff09; 有效 IP 地址 正好由四个整数&#xff08;每个整数位于 0 到 255 之间组成&#xff0c;且不能含有前导 0&#xff09;&#xff0c;整数之间用 . 分隔。 例如&#xff1a;"0.1.2.201" 和 "192.168.1.1…

CS 2520nonono

CS 2520nonono WeChat&#xff1a;yj4399_​​​​​ Sina Visitor System High-level●3 Congestion Control Algorithms:○TCP Reno:■additive increase, multiplicative decrease function to adjust window size for every RTTuntil a packet loss is detected○TCP CUBI…

用java实现拼图小游戏

1、了解拼图游戏基本功能&#xff1a; 拼图游戏内容由若干小图像块组成的&#xff0c;通过鼠标点击图像块上下左右移动&#xff0c;完成图像的拼凑。 2、拼图游戏交互界面设计与开发&#xff1a; 通过创建窗体类、菜单、中间面板和左右面板完成设计拼图的交互界面 &#xff…

分享从零开始学习网络设备配置--任务4.3 使用动态路由RIPng实现网络连通

任务描述 某公司使用IPv6技术搭建企业网络&#xff0c;由于静态路由需要管理员手工配置&#xff0c;在网络拓扑发生变化时&#xff0c;也不会自动生成新的路由&#xff0c;因此采用IPv6动态路由协议RIPng实现网络连通&#xff0c;实现任意两个节点之间的通信&#xff0c;并降低…

基于SpringBoot学生读书笔记共享

摘 要 本论文主要论述了如何使用JAVA语言开发一个读书笔记共享平台 &#xff0c;本系统将严格按照软件开发流程进行各个阶段的工作&#xff0c;采用B/S架构&#xff0c;面向对象编程思想进行项目开发。在引言中&#xff0c;作者将论述读书笔记共享平台的当前背景以及系统开发的…

第16关 革新云计算:如何利用弹性容器与托管K8S实现极速服务POD扩缩容

------> 课程视频同步分享在今日头条和B站 天下武功&#xff0c;唯快不破&#xff01; 大家好&#xff0c;我是博哥爱运维。这节课给大家讲下云平台的弹性容器实例怎么结合其托管K8S&#xff0c;使用混合服务架构&#xff0c;带来极致扩缩容快感。 下面是全球主流云平台弹…

对抗产品团队中的认知偏误:给产品经理的专家建议

今天的产品经理面临着独特的挑战。他们不仅需要设计和构建创新功能&#xff0c;还必须了解这些功能将如何为客户带来价值并推进关键业务目标。如果不加以控制&#xff0c;认知偏差可能会导致您构建的内容与客户想要的内容或业务需求之间不一致。本文将详细阐述产品经理可以避免…

下载MySQL JDBC驱动的方法

说明 java代码通过JDBC访问MySQL数据库&#xff0c;需要MySQL JDBC驱动。 例如&#xff0c;下面这段代码&#xff0c;因为找不到JDBC驱动&#xff0c;所以执行会报异常&#xff1a; package com.thb;public class JDBCDemo {public static void main(String[] args) throws …

网络基础_1

目录 网络基础 协议 协议分层 OSI七层模型 网络传输的基本流程 数据包的封装和分用 IP地址和MAC地址 网络基础 网络就是不同的计算机之间可以进行通信&#xff0c;前面我们学了同一台计算机之间通信&#xff0c;其中有进程间通信&#xff0c;前面学过的有管道&#xff…