R语言：SVD分解求解线性方程组AX=b

2024-04-11 22:01:20

R语言：SVD分解求解线性方程组AX=b

当AAA是方阵时，可以直接用内置函数解，当A" role="presentation" style="position: relative;">AAA不是方阵时，只能求得最小二乘解。

函数svd的用法

svd分解X，使用函数svd，返回一个列表T，顺序是d, u, v。

X=UDV′X=UDV′

X=UDV'

T <- svd(X)
U <- T$u
D <- T$d
V <- T$v

这里T是list，注意这里的U和V是矩阵，D是向量，想要恢复原矩阵X，需要：

Y <- U %*% diag(D) %*% t(V)

求解方法

问题：求xxx使得(最小二乘解)：

线性方程组形式

Ax=b" role="presentation" style="text-align: center; position: relative;">Ax=bAx=b

Ax=b
对AAA进行SVD分解：

A=UDV′" role="presentation" style="text-align: center; position: relative;">A=UDV′A=UDV′

A=UDV'
AAA的广义逆矩阵A+" role="presentation" style="position: relative;">A+A+A^+：

A+=VD+UA+=VD+U

A^+=VD^+U
其中 D+D+D^+是 DDD的广义逆：
广义逆的定义是：

D+=(D‘D)−1D′" role="presentation" style="text-align: center; position: relative;">D+=(D‘D)−1D′D+=(D‘D)−1D′

D^+=(D‘D)^{-1}D’
数值分析书上讲，是非零元的逆的转置。但此时计算时出现问题：

> A <- matrix(rep(1,18),nrow = 3)
#### A的秩为1
> T <- svd(A)
#### 对A进行svd分解
> A2 <- T$u %*% diag(T$d) %*% t(T$v)
#### 利用svd分解恢复A
> diag(T$d)
# [,1]         [,2]         [,3]
# [1,] 4.242641 0.000000e+00 0.000000e+00
# [2,] 0.000000 8.107923e-16 0.000000e+00
# [3,] 0.000000 0.000000e+00 6.972611e-32

这里就出现问题了，虽然A的秩为1，但是计算得到的奇异值，由于精度原因，后两个为很接近于0的数。

所以要利用svd算A的广义逆，需要忽略后两项。取diag(T$d)的前r项，这个r就是A的rank。求A的rank可以用qr分解

r <- qr(A)$rank

在数值计算中，太小的奇异值会被忽略，这就成为低秩分解。

还有一种方法，直接设置阈值θθ\theta，小于θθ\theta的都记为0，但是这个θθ\theta的取法我还没有查到文献记载。

那么：

x=A+bx=A+b

x =A^+b

矩阵形式

AX=BAX=B

AX=B
则

X=A+BX=A+B

X=A^+B

系数矩阵为右乘

XA=BXA=B

XA=B
对上式求转置即可：

A′X′=B′A′X′=B′

A'X'=B'
对 A′A′A'进行SVD分解，按照上面计算得到 X′X′X'，再转置得到 XX<script type="math/tex" id="MathJax-Element-2830">X</script>。

R语言：SVD分解求解线性方程组AX=b相关推荐

共轭梯度法求解线性方程组Ax=b（附代码）
共轭梯度法求解线性方程组:Ax=b 数值分析老师给的作业,写出来代码平时分满分.简单勿喷,hhhh 原理 CG求解线性方程组的原理大家翻一翻数值分析的课本即可,此处不哔哔直接上matlab代码代码 ...
python QR分解求解线性方程组和矩阵本征值和本征向量
下面的代码提供了两个函数 solve_linear_equ, 利用QR分解求解线性方程组,输入是一个二维的非奇异的系数方阵和一个常数array,输出是该线性方程组的解 eigen, 输入是一个实方阵( ...
课程设计—C++实现高斯消元法求解线性方程组Ax=b（附源码）
(一)需求和规格说明输入是 N(N<256) 元线性方程组 Ax=B, 输出是方程组的解,也可能无解或有多组解.可以用高斯消去法求解,也可以采用其它方法. 要求给出线性方程组的矩阵,能够输出线 ...
基于R语言SVD的图像压缩方法
图片的基本构成假设我们有一个灰度图像(128×128,即128×128像素).我们可以使用矩阵来表示这个图像.如果我们有彩色图像,用R中的readJEPG函数读取图片就可以发现它是由一个三维数组构成 ...
fortran基于svd分解求解广义逆矩阵
fortran的MKl函数库中好像没有求解广义逆矩阵的sub 本篇文章给出基于svd分解求广义逆矩阵的示例代码 program pinvuse lapack95implicit noneinteger ...
c语言如何计算出迭代次数,计算方法——C语言实现——迭代法求解线性方程组...
最近在上计算方法这门课,要求是用MATLAB做练习题,但是我觉得C语言也很棒棒啊~ 题目: 和直接法不同,迭代法是一种逐次逼近的方法,将复杂问题简单化,求比较大的方程组时一般都不会用直接法.迭代法有好 ...
数值分析——LU分解求解线性方程组的Python实现
import numpy as np import math A=np.array([[1,2,3],[2,5,2],[3,1,5]]) # np.mat创建矩阵,np.arry创建数组 b=np.a ...
matlab ax=b求x,Matlab求解线性方程组Ax=b的几种常见方法Matlab求解线性方程组Ax=b的几种常见方法...
例如方程组: 法1:左除法 >> A=[3 1 -1;1 2 4;-1 4 5];b=[3.6;2.1;-1.4]; >> x=A\b x = 1.4818 -0.4606 0 ...
高斯消元法解矩阵方程c语言,c++高斯消元法求解线性方程组
//高斯消元法求解方程组 #include #include using namespace std; #define MaxNum 10 int array[MaxNum][MaxNum] = { ...
解方程AX=b与矩阵分解：奇异值分解（SVD分解）特征值分解 QR分解三角分解 LLT分解
文章目录 1. 前言 2. LU三角分解 3. Cholesky分解 - LDLT分解 4. Cholesky分解 - LLT分解 5. QR分解 6. 奇异值分解 7. 特征值分解 1. 前言本博 ...

最新文章

热门文章