MATLAB中qr函数用法

目录

语法

说明

示例

Q-Less QR 分解

矩阵的完整 QR 分解

置换 QR 分解

用精简 QR 因子求解线性系统

求解稀疏线性系统

求解矩形稀疏线性系统

提示


        qr函数的功能是对矩阵进行QR 分解。

语法

R = qr(A)
[Q,R] = qr(A)
[Q,R,P] = qr(A)
[___] = qr(A,"econ")
[Q,R,P] = qr(A,outputForm)
[___] = qr(A,0)
[C,R] = qr(S,B)
[C,R,P] = qr(S,B)
[___] = qr(S,B,"econ")
[C,R,P] = qr(S,B,outputForm)
[___] = qr(S,B,0)

说明

R = qr(A) 返回 QR 分解 A = Q*R 的上三角 R 因子。

[Q,R] = qr(A) 对 m×n 矩阵 A 执行 QR 分解,满足 A = Q*R。因子 R 是 m×n 上三角矩阵,因子 Q 是 m×m 正交矩阵。

[Q,R,P] = qr(A) 还返回一个置换矩阵 P,满足 A*P = Q*R。如果 A 为满矩阵,将选择置换矩阵,使得 abs(diag(R)) 递减。

[___] = qr(A,"econ") 使用上述任意输出参数组合进行精简分解。输出的大小取决于 m×n 矩阵 A 的大小:

        如果 m > n,则 qr 仅计算 Q 的前 n 列和 R 的前 n 行。

        如果 m <= n,则精简分解与常规分解相同。

[Q,R,P] = qr(A,outputForm) 指定置换信息 P 是以矩阵还是向量形式返回。例如,如果 outputForm 是 "vector",则 A(:,P) = Q*R。outputForm 的默认值是 "matrix",满足 A*P = Q*R。

[___] = qr(A,0) 等效于 qr(A,"econ","vector")。不建议使用此语法。改用 "econ" 选项。

[C,R] = qr(S,B) 计算 C = Q'*B 和上三角因子 R。您可以使用 C 和 R 计算稀疏线性系统 S*X = B 和 X = R\C 的最小二乘解。

[C,R,P] = qr(S,B) 还返回置换矩阵 P,选择该矩阵是为了减少 R 中的填充。您可以使用 C、R 和 P 计算稀疏线性系统 S*X = B 和 X = P*(R\C) 的最小二乘解。

[___] = qr(S,B,"econ") 使用上述任意输出参数组合进行精简分解。输出的大小取决于 m×n 稀疏矩阵 S 的大小:

        如果 m > n,则 qr 仅计算 C 和 R 的前 n 行。

        如果 m <= n,则精简分解与常规分解相同。

[C,R,P] = qr(S,B,outputForm) 指定置换信息 P 是以矩阵还是向量形式返回。例如,如果 outputForm 是 "vector",则 S*X = B 的最小二乘解是 X(P,:) = R\C。outputForm 的默认值为 "matrix",此时 S*X = B 的最小二乘解是 X = P*(R\C)。

[___] = qr(S,B,0) 等效于 qr(S,B,"econ","vector")。不建议使用此语法。改用 "econ" 选项。

示例

Q-Less QR 分解

        求 5×5 幻方矩阵的 QR 分解。指定一个输出参数以只返回上三角因子。

A = magic(5);
R = qr(A)
R = 5×5

  -32.4808  -26.6311  -21.3973  -23.7063  -25.8615
         0   19.8943   12.3234    1.9439    4.0856
         0         0  -24.3985  -11.6316   -3.7415
         0         0         0  -20.0982   -9.9739
         0         0         0         0  -16.0005

矩阵的完整 QR 分解

        通过指定两个输出参数来计算幻方测试矩阵的完整 QR 分解。

A = magic(5);
[Q,R] = qr(A)
Q = 5×5

   -0.5234    0.5058    0.6735   -0.1215   -0.0441
   -0.7081   -0.6966   -0.0177    0.0815   -0.0800
   -0.1231    0.1367   -0.3558   -0.6307   -0.6646
   -0.3079    0.1911   -0.4122   -0.4247    0.7200
   -0.3387    0.4514   -0.4996    0.6328   -0.1774

R = 5×5

  -32.4808  -26.6311  -21.3973  -23.7063  -25.8615
         0   19.8943   12.3234    1.9439    4.0856
         0         0  -24.3985  -11.6316   -3.7415
         0         0         0  -20.0982   -9.9739
         0         0         0         0  -16.0005

        在计算机精度范围内验证 A=QR。

norm(A-Q*R)
ans = 1.2642e-14

置换 QR 分解

        指定三个输出参数以返回一个置换矩阵或向量,该置换矩阵或向量可减少 QR 分解的 R 因子中的填充。

        计算 west0479 稀疏矩阵的 QR 分解。指定三个输出以返回满足 AP=QR 的置换矩阵。

load west0479
A = west0479;
[Q,R,P] = qr(A);

        在计算机精度范围内验证置换矩阵 P 满足 A*P = Q*R。

norm(A*P-Q*R,"fro")
ans = 3.5411e-10

        现在指定 "vector" 选项以将 p 以置换向量形式返回。

[Q,R,p] = qr(A,"vector");

        在计算机精度范围内验证置换向量 p 满足 A(:,p) = Q*R。

norm(A(:,p) - Q*R,"fro")
ans = 3.5411e-10

        验证对于稀疏输入,在分解中使用置换矩阵或置换向量所得的R因子的非零项比使用非置换分解的少。

[Q1,R1] = qr(A);
spy(R1)

如图所示:

spy(R)

如图所示:

        结果表明置换分解产生的 R 因子的非零值明显减少。

用精简 QR 因子求解线性系统

        使用系数矩阵的精简 QR 分解来求解线性系统 Ax=b。使用 magic(10) 的前五列创建一个 10×5 系数矩阵。对于线性方程 Ax=b 的右侧,使用矩阵的行总和。在这种设置下,方程 x 的解应为由 1 组成的向量。

A = magic(10);
A = A(:,1:5)
A = 10×5

    92    99     1     8    15
    98    80     7    14    16
     4    81    88    20    22
    85    87    19    21     3
    86    93    25     2     9
    17    24    76    83    90
    23     5    82    89    91
    79     6    13    95    97
    10    12    94    96    78
    11    18   100    77    84

b = sum(A,2)
b = 10×1

   215
   215
   215
   215
   215
   290
   290
   290
   290
   290

        计算 A 的精简 QR 分解。然后用 x(p,:) = R\(Q\b) 求解线性系统 QRx=b。由于 Q 是正交矩阵,此方程与 x(p,:) = R\(Q'*b) 相同。

[Q,R,p] = qr(A,"econ","vector")
Q = 10×5

   -0.0050   -0.4775   -0.0504    0.5193    0.0399
   -0.0349   -0.5001   -0.0990   -0.1954   -0.2006
   -0.4384    0.1059   -0.4660    0.4464    0.0628
   -0.0947   -0.4151   -0.2923   -0.2542    0.5274
   -0.1246   -0.4117   -0.2812   -0.1326   -0.4130
   -0.3787    0.0209    0.2702    0.4697    0.0390
   -0.4085   -0.0017    0.2217   -0.2450   -0.2015
   -0.0648   -0.3925    0.6939    0.0669    0.1225
   -0.4683    0.0833    0.0283   -0.3038    0.5265
   -0.4982    0.0867    0.0394   -0.1822   -0.4138

R = 5×5

 -200.7112  -55.5026 -167.6040  -84.7237 -168.7997
         0 -192.1053  -40.3557 -152.4040  -39.2814
         0         0  101.3180  -89.4254   96.0172
         0         0         0   41.0248  -14.9083
         0         0         0         0   24.6386

p = 1×5

     3     1     5     2     4

x(p,:) = R\(Q\b)
x = 5×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000

        生成 R 的对角线的半对数图,以确认置换分解产生的 R 因子的 abs(diag(R)) 是递减的。将 A 的奇异值绘制在同一个图中进行比较。实际上,R 的对角线值与 A 的奇异值具有类似的行为。因此,可以使用 R 的对角线值来衡量矩阵 A 的奇异程度。

semilogy(abs(diag(R)),"-o")
hold on
semilogy(svd(A),"r-o")
legend("Diagonal of R","Singular Values of A")

如图所示:

求解稀疏线性系统

        求解稀疏线性系统,并使用结果查看向量 b 有多少位于 S 的列空间中。创建一个随机 500×20 稀疏矩阵(密度为 10%)和一个由 1 组成的向量。使用 qr 将矩阵分解为因子 R 和 C = Q'*b

S = sprand(500,20,0.1);
b = ones(500,1);
[C,R] = qr(S,b,"econ");

        使用这些结果求解 Sx=b,解为 x = R\C。

x = R\C;

        以恒等式 为例。除以 b 的范数,会得到一个新的恒等式,表明 b 有多少位于 S 的列空间中:

        第一个项说明 b 有多少不在 S 的列空间中,而第二个项说明 b 有多少在 S 的列空间中。

t1 = norm(S*x-b)^2/norm(b)^2
t1 = 0.4000
t2 = norm(C)^2/norm(b)^2
t2 = 0.6000

求解矩形稀疏线性系统

        使用 qr 求解矩阵方程 Sx=B,其中 S 为矩形稀疏系数矩阵。

        加载 west0479 稀疏矩阵,并将前 200 列用作线性系统中的矩形系数矩阵。对于方程的右侧,使用 S 的行总和。在这种设置下,Sx=B 的解是由 1 组成的向量。

load west0479
S = west0479(:,1:200);
B = sum(S,2);

        使用 qr 以及两个输入和三个输出求解 Sx=B。线性系统的解是 x = P*(R\C)。

[C,R,P] = qr(S,B);
x = P*(R\C);

在计算机精度范围内验证 Sx−B=0。

norm(S*x-B)
ans = 9.1703e-11

        注意:要计算上三角因子 R 和置换矩阵 P,但避免计算正交矩阵 Q(这通常是对 qr 的调用中计算量最大的部分),可以将 B 指定为空矩阵:

emptyB = zeros(size(S,1),0);
[~,R,P] = qr(S,emptyB);

提示

  • ​要求解涉及相同系数矩阵的多个线性系统,请使用 decomposition 对象。

  • 对于语法 [C,R] = qr(S,B),仅当 S 不具有低秩时,X = R\C 的值才是 S*X = B 的最小二乘解。

  • 17
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
### 回答1: polyfit函数MATLAB的一个用于多项式拟合的函数,它在MATLAB的polyfit文档有详细的说明。但是,polyfit的具体实现是在MATLAB的底层通过C语言编写的。由于MATLAB的源代码是闭源的,因此无法直接获取polyfit函数的C语言源代码。 polyfit函数算法通常使用最小二乘法来拟合多项式曲线。它通过计算数据点与拟合曲线之间的误差,并调整多项式系数,使误差最小化。MATLAB的polyfit函数使用的算法可能是基于一些经典的数值方法,如QR分解或LU分解。 如果你需要在C语言实现多项式拟合的功能,可以参考polyfit函数算法步骤,并使用C语言编写相应的代码。您可以根据多项式拟合的算法自行实现,如使用最小二乘法或其他数值方法,通过计算误差最小化来调整多项式系数。 这里是一个简单的示例,使用C语言实现2次多项式拟合的算法步骤: 1. 输入原始数据点的坐标(x, y)。 2. 定义拟合多项式的阶数n为2。 3. 创建一个矩阵A和一个向量B,用于存储方程组Ax=B的系数和常数项。 4. 遍历所有的数据点,分别计算矩阵A和向量B的元素。 - A的第(i, j)个元素是x的i次方的和,其j表示多项式的次数。 - B的第(i)个元素是y与x的i次方的乘积的和。 5. 使用数值方法(如高斯消元法或QR分解)求解方程组Ax=B,得到多项式系数。 6. 输出多项式的系数,即将多项式曲线用一组Cn表示。 需要注意的是,以上是一个简化的示例,实际的C语言实现可能需要更多的代码和复杂的数值计算方法。如果你需要更详细的C语言源代码实现,可以参考相关的数值计算库,如GSL(GNU Scientific Library)或使用其他开源的数值计算库来实现多项式拟合。 ### 回答2: polyfit函数MATLAB用于多项式拟合的函数,用于通过最小二乘法拟合数据点到一个多项式曲线。polyfit函数的C语言源代码如下: ```c #include <stdio.h> #include <stdlib.h> void polyfit(double *x, double *y, int n, double *coefficients, int order) { int i, j, k; double *X = (double *) malloc((2 * order + 1) * sizeof(double)); double *Y = (double *) malloc((order + 1) * sizeof(double)); double *B = (double *) malloc((order + 1) * sizeof(double)); double *A = (double *) malloc((order + 1) * (order + 1) * sizeof(double)); for (i = 0; i < 2 * order + 1; i++) { X[i] = 0.0; for (j = 0; j < n; j++) { X[i] += pow(x[j], i); } } for (i = 0; i <= order; i++) { Y[i] = 0.0; for (j = 0; j < n; j++) { Y[i] += pow(x[j], i) * y[j]; } } for (i = 0; i <= order; i++) { for (j = 0; j <= order; j++) { A[i * (order + 1) + j] = X[i + j]; } } for (i = 0; i <= order; i++) { B[i] = Y[i]; } for (k = 0; k <= order; k++) { for (i = k + 1; i <= order; i++) { double factor = A[i * (order + 1) + k] / A[k * (order + 1) + k]; B[i] -= factor * B[k]; for (j = k; j <= order; j++) { A[i * (order + 1) + j] -= factor * A[k * (order + 1) + j]; } } } coefficients[order] = B[order] / A[order * (order + 1) + order]; for (i = order - 1; i >= 0; i--) { double sum = B[i]; for (j = i + 1; j <= order; j++) { sum -= A[i * (order + 1) + j] * coefficients[j]; } coefficients[i] = sum / A[i * (order + 1) + i]; } free(X); free(Y); free(B); free(A); } int main() { double x[] = {1, 2, 3, 4, 5}; double y[] = {2, 3, 5, 8, 10}; int n = 5; int order = 2; double *coefficients = (double *) malloc((order + 1) * sizeof(double)); polyfit(x, y, n, coefficients, order); for (int i = 0; i <= order; i++) { printf("Coefficient %d: %.2f\n", i, coefficients[i]); } free(coefficients); return 0; } ``` 这是一个简单的多项式拟合的例子,输入的数据点为x和y,n为数据点个数,order为拟合多项式的阶数。在主函数调用polyfit函数进行拟合,拟合结果存储在coefficients数组,然后打印出每个系数的值。 需要注意的是,此为简化版本的多项式拟合代码,实际情况可能还需要添加其他处理和优化策略,以适应更加复杂和实际的数据拟合需求。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值