分享

矩阵

 圣树荫 2013-07-05

矩阵的基本运算

§3.1 加和减

如矩阵AB的维数相同,则A+BA-B表示矩阵AB的和与差.对应位置相加减

A=                                B=

1     2     3                   1     4     7

4     5     6                   2     5     8

7     8     0                   3     6     0

C =A+B返回:

C =

     2     6    10

     6    10    14

    10    14     0

如果运算对象是个标量(即1×1矩阵),可和其它矩阵进行加减运算.例如:

x=    -1                  y=x-1=    -2

0                                 -1

2                                 1                  

§3.2矩阵乘法

§3.2.1 矩阵的普通乘法

矩阵乘法用“ * ”符号表示,当A矩阵列数与B矩阵的行数相等时,二者可以进行乘法运算,否则是错误的.计算方法和线性代数中所介绍的完全相同. i行与j列个数字依次相乘之和为结果的“Aij”


§3.3 矩阵除法

Matlab中有两种矩阵除法符号:“\”即左除和“/”即右除.如果A矩阵是非奇异方阵,则A\BA的逆矩阵乘B,即inv(A)*B;而B/ABA的逆矩阵,即B*inv(A).具体计算时可不用逆矩阵而直接计算.通常:x=A\B就是A*x=B的解;

x=B/A就是x*A=B的解.

BA矩阵行数相等可进行左除.如果A是方阵,用高斯消元法分解因数.解方程:A*x(:, j)=B(:, j),式中的(:, j)表示B矩阵的第j列,返回的结果x具有与B矩阵相同的阶数,如果A是奇异矩阵将给出警告信息.

如果A矩阵不是方阵,可由以列为基准的Householder正交分解法分解,这种分解法可以解决在最小二乘法中的欠定方程或超定方程,结果是m×nx矩阵.mA矩阵的列数,nB矩阵的列数.每个矩阵的列向量最多有k个非零元素,k A的有效秩.

右除B/A可由B/A=(A'\B')'左除来实现.

§3.4矩阵乘方

A^P意思是AP次方.如果A是一个方阵,P是一个大于1的整数,则A^P表示AP次幂,即A自乘P次.如果P不是整数,计算涉及到特征值和特征向量的问题,如已经求得:[V,D]=eig(A),则:

A^P=V*D.^P/V(注:这里的.^表示数组乘方,或点乘方,参见后面的有关介绍)

如果B是方阵, a是标量,a^B就是一个按特征值与特征向量的升幂排列的B次方程阵. 如果aB都是矩阵,则a^B是错误的

§3.5 矩阵的超越函数

Matlab中解释exp(A)sqrt(A)时曾涉及到级数运算,此运算定义在A的单个元素上. Matlab可以计算矩阵的超越函数,如矩阵指数、矩阵对数等.

一个超越函数可以作为矩阵函数来解释,例如将“m”加在函数名的后边而成expm(A)sqrtm(A),当Matlab运行时,有下列三种函数定义:

expm           矩阵指数

logm        矩阵对数

sqrtm           矩阵开方

所列各项可以加在多种m文件中或使用funm.请见应用库中sqrtm.m1ogm.mfunm.m文件和命令手册.

§3.6数组运算

数组运算由线性代数的矩阵运算符“*”、“/”、“\”、“^”前加一点来表示,即为“.*”、“./”、“.\”、“.^”.注意没有“.+”、“.-”运算

§3.6.1数组的加和减

对于数组的加和减运算与矩阵运算相同,所以“+”、“-”既可被矩阵接受又可被数组接受.

§3.6.2数组的乘和除

数组的乘用符号.*表示,如果AB矩阵具有相同阶数,则A.*B表示AB单个元素之间的对应相乘.例如 x=[1  2  3]; y=[ 4  5  6];

计算z=x.*y

结果z=4  10  18

数组的左除(.\)与数组的右除(./),由读者自行举例加以体会.

§3.6.3 数组乘方

数组乘方用符号.^表示.

例如:键入:

x=[ 1  2  3]

y=[ 4  5  6]

z=x.^y=[1^4  2^5  3^6]=[1  32  729]

(1) 如指数是个标量,例如x.^2x同上,则:

z=x.^2=[1^2  2^2  3^2]=[ 1  4  9]

(2) 如底是标量,例如2 .^[x y] xy同上,则:

z=2 .^[x y]=[2^1  2^2  2^3  2^4  2^5  2^6]=[2  4  8  16  32  64]

从此例可以看出Matlab算法的微妙特性,虽然看上去与其它乘方没什么不同,但在2和“.”之间的空格很重要,如果不这样做,解释程序会把“.”看成是2的小数点. Matlab看到符号“^”时,就会当做矩阵的幂来运算,这种情况就会出错,因为指数矩阵不是方阵.

§3.7 矩阵函数

Matlab的数学能力大部分是从它的矩阵函数派生出来的,其中一部分装入Matlab本身处理中,它从外部的Matlab建立的M文件库中得到,还有一些由个别的用户为其自己的特殊的用途加进去的.其它功能函数在求助程序或命令手册中都可找到.手册中备有为Matlab提供数学基础的LINPACKEISPACK软件包,提供了下面四种情况的分解函数或变换函数:

(1)三角分解;(2)正交变换;(3) 特征值变换;(4)奇异值分解.

§3.7.1三角分解

最基本的分解为“LU”分解,矩阵分解为两个基本三角矩阵形成的方阵,三角矩阵有上三角矩阵和下三角矩阵.计算算法用高斯变量消去法.

lu函数中可以得到分解出的上三角与下三角矩阵,函数inv得到矩阵的逆矩阵,det得到矩阵的行列式.解线性方程组的结果由方阵的“\”和“/”矩阵除法来得到.

例如:

A=[ 1      2     3

4      5     6

7      8     0]

LU分解,用Matlab的多重赋值语句

[L,U]=lu(A)

得出

L =

0.1429

1.0000

0

0.5714

0.5000

1.0000

1.0000

0

0

U =

7.0000

8.0000

0

0

0.8571

3.0000

0

0

4.5000

注:L是下三角矩阵的置换,U是上三角矩阵的正交变换,分解作如下运算,检测计算结果只需计算L*U即可.

求逆由下式给出: x=inv(A)

x =

-1.7778

0.8889

-0.1111

1.5556

-0.7778

0.2222

-0.1111

0.2222

-0.1111

LU分解得到的行列式的值是精确的,d=det(U)*det(L)的值可由下式给出:

d=det(A)

d =

27

直接由三角分解计算行列式:d=det(L)*det(U)

d =

27.0000

为什么两种d的显示格式不一样呢? Matlabdet(A)运算时,所有A的元素都是整数,所以结果为整数.但是用LU分解计算d时,LU的元素是实数,所以Matlab产生的d也是实数.

例如:线性联立方程取 b=[ 1

                                           3

                                           5]

Ax=b方程,用Matlab矩阵除得到

x=A\b

结果x=

0.3333

0.3333

0.0000

由于A=L*U,所以x也可以有以下两个式子计算:y=L\bx=U\y.得到相同的x值,中间值y为:

y =

5.0000

0.2857

0.0000

Matlab中与此相关的函数还有rcondcholrref.其基本算法与LU分解密切相关.chol函数对正定矩阵进行Cholesky分解,产生一个上三角矩阵,以使R'*R=Xrref用具有部分主元的高斯-约当消去法产生矩阵A的化简梯形形式.虽然计算量很少,但它是很有趣的理论线性代数.为了教学的要求,也包括在Matlab中.

§3.7.2正交变换

QR”分解用于矩阵的正交-三角分解.它将矩阵分解为实正交矩阵或复酉矩阵与上三角矩阵的积,对方阵和长方阵都很有用.

例如A=[  1     2     3

4     5     6

7     8     9

10    11    12]

是一个降秩矩阵,中间列是其它二列的平均,我们对它进行QR分解:

[Q,R]=qr(A)

Q =

-0.0776

-0.8331

0.5444

0.0605

-0.3105

-0.4512

-0.7709

0.3251

-0.5433

-0.0694

-0.0913

-0.8317

-0.7762

0.3124

0.3178

0.4461

R =

-12.8841

-14.5916

-16.2992

0

-1.0413

-2.0826

0

0

0.0000

0

0

0

可以验证Q*R就是原来的A矩阵.由R的下三角都给出0,并且R(3,3)=0.0000,说明矩阵R与原来矩阵A都不是满秩的.

下面尝试利用QR分解来求超定和降秩的线性方程组的解.

例如:

b=[ 1

3

5

7]

讨论线性方程组Ax=b,我们可以知道方程组是超定的,采用最小二乘法的最好结果是计算x=A\b

结果为:

Warning: Rank deficient, rank = 2  tol =   1.4594e-014

x =

    0.5000

         0

    0.1667

我们得到了缺秩的警告.用QR分解法计算此方程组分二个步骤:

y=Q'*b

x=R\y

求出的y值为

y =

-9.1586

-0.3471

0.0000

0.0000

x的结果为

Warning: Rank deficient, rank = 2  tol =   1.4594e-014

x =

    0.5000

         0

    0.1667

A*x来验证计算结果,我们会发现在允许的误差范围内结果等于b.这告诉我们虽然联立方程Ax=b是超定和降秩的,但两种求解方法的结果是一致的.显然x向量的解有无穷多个,而“QR”分解仅仅找出了其中之一.

§3.7.3奇异值分解

Matlab中三重赋值语句

[U,S,V]=svd(A)

在奇异值分解中产生三个因数:

A=U*S*V '

U矩阵和V矩阵是正交矩阵,S矩阵是对角矩阵,svd(A)函数恰好返回S的对角元素,而且就是A的奇异值(其定义为:矩阵A'*A的特征值的算术平方根).注意到A矩阵可以不是方的矩阵.

奇异值分解可被其它几种函数使用,包括广义逆矩阵pinv(A)、秩rank(A)、欧几里德矩阵范数norm(A,2)和条件数cond(A)

§3.7.4 特征值分解

如果An×n矩阵,若l满足Ax=lx,则称lA的特征值,x为相应的特征向量.

函数eig(A)返回特征值列向量,如果A是实对称的,特征值为实数.特征值也可能为复数,例如:

A=[   0     1

-1     0]

eig(A)

产生结果

ans =

0 + 1.0000i

0 - 1.0000i

如果还要求求出特征向量,则可以用eig(A)函数的第二个返回值得到:

[x,D]=eig(A)

D的对角元素是特征值.x的列是相应的特征向量,以使A*x=x*D

计算特征值的中间结果有两种形式:

Hessenberg形式为hess(A)Schur形式为schur(A)

schur形式用来计算矩阵的超越函数,诸如sqrtm(A)logm(A)

如果AB是方阵,函数eig(A,B)返回一个包含一般特征值的向量来解方程

Ax=lBx

双赋值获得特征向量

[X,D]=eig(A,B)

产生特征值为对角矩阵D.满秩矩阵X的列相应于特征向量,使A*X=B*X*D,中间结果由qz(A,B)提供

§3.7.5

Matlab计算矩阵A的秩的函数为rank(A),与秩的计算相关的函数还有:rref(A)orth(A)null(A)和广义逆矩阵pinv(A)等.

利用rref(A)A的秩为非0行的个数.rref方法是几个定秩算法中最快的一个,但结果上并不可靠和完善.pinv(A)是基于奇异值的算法.该算法消耗时间多,但比较可靠.其它函数的详细用法可利用Help求助.

                                                                                                                                   


    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。
    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多