当前位置：文档库 › MatCalc----矩阵计算

MatCalc----矩阵计算

Notes on Matrix Calculus

Paul L.Fackler?

North Carolina State University

September27,2005

Matrix calculus is concerned with rules for operating on functions of matrices.For example,suppose that an m×n matrix X is mapped into a p×q matrix Y.We are interested in obtaining expressions for derivatives such as

?Y ij

?X kl

for all i,j and k,l.The main di?culty here is keeping track of where things are put.There is no reason to use subscripts;it is far better instead to use a system for ordering the results using matrix operations.

Matrix calculus makes heavy use of the vec operator and Kronecker products.The vec operator vectorizes a matrix by stacking its columns(it is convention that column rather than row stacking is used).For example, vectorizing the matrix

12 34 56

?Paul L.Fackler is an Associate Professor in the Department of Agricultural and Re-source Economics at North Carolina State University.These notes are copyrighted mate-rial.They may be freely copied for individual use but should be appropriated referenced in published work.

Mail:Department of Agricultural and Resource Economics

NCSU,Box8109

Raleigh NC,27695,USA

e-mail:paul fackler@https://www.wendangku.net/doc/171169668.html,

Web-site:https://www.wendangku.net/doc/171169668.html,/～pfackler/

c 2005,Paul L.Fackler

produces

The Kronecker product of two matrices,A and B,where A is m×n and B is p×q,is de?ned as

A?B=

A11B A12B...A1n B

A21B A22B...A2n B

............

A m1

B A m2B...A mn B

which is an mp×nq matrix.There is an important relationship between the Kronecker product and the vec operator:

vec(AXB)=(B ?A)vec(X).

This relationship is extremely useful in deriving matrix calculus results.

Another matrix operator that will prove useful is one related to the vec operator.De?ne the matrix T m,n as the matrix that transforms vec(A)into vec(A ):

T m,n vec(A)=vec(A ).

Note the size of this matrix is mn×mn.T m,n has a number of special properties.The?rst is clear from its de?nition;if T m,n is applied to the vec of an m×n matrix and then T n,m applied to the result,the original vectorized matrix results:

T n,m T m,n vec(A)=vec(A).

Thus

T n,m T m,n=I mn.

The fact that

T n,m=T?1

m,n

follows directly.Perhaps less obvious is that

T m,n=T n,m

(also combining these results means that T m,n is an orthogonal matrix).

The matrix operator T m,n is a permutation matrix,i.e.,it is composed of0s and1s,with a single1on each row and column.When premultiplying another matrix,it simply rearranges the ordering of rows of that matrix (postmultiplying by T m,n rearranges columns).

The transpose matrix is also related to the Kronecker product.With A and B de?ned as above,

B?A=T p,m(A?B)T n,q.

This can be shown by introducing an arbitrary n×q matrix C:

T p,m(A?B)T n,q vec(C)=T p,m(A?B)vec(C )

=T p,m vec(BC A )

=vec(ACB )

=(B?A)vec(C).

This implies that((B?A)?T p,m(A?B)T n,q)vec(C)=0.Because C is arbitrary,the desired result must hold.

An immediate corollary to the above result is that

(A?B)T n,q=T m,p(B?A).

It is also useful to note that T1,m=T m,1=I m.Thus,if A is1×n then (A?B)T n,q=(B?A).When working with derivatives of scalars this can result in considerable simpli?cation.

Turning now to calculus,de?ne the derivative of a function mapping n→ m as the m×n matrix of partial derivatives:

[Df]ij=?f i(x)?x j

For example,the simplest derivative is

dAx

=A.

Using this de?nition,the usual rules for manipulating derivatives apply naturally if one respects the rules of matrix conformability.The summation rule is obvious:

D[αf(x)+βg(x)]=αDf(x)+βDg(x),

whereαandβare scalars.The chain rule involves matrix multiplication, which requires conformability.Given two functions f: n→ m and g: p→ n,the derivative of the composite function is

D[f(g(x))]=f (g(x))g (x).

Notice that this satis?es matrix multiplication conformability,whereas the expression g (x)f (g(x))attempts to postmultipy an n×p matrix by an m×n matrix.To de?ne a product rule,consider the expression f(x) g(x), where f,g: n→ m.The derivative is the1×n vector given by

D[f(x) g(x)]=g(x) f (x)+f(x) g (x).

Notice that no other way of multiplying g by f and f by g would ensure conformability.A more general version of the product rule is de?ned below.

The product rule leads to a useful result about quadratic functions:

dx Ax

=x A+x A =x (A+A ).

When A is symmetric this has the very natural form dx Ax/dx=2x A.

These rules de?ne derivatives for vectors.De?ning derivatives of matri-ces with respect to matrices is accomplished by vectorizing the matrices,so dA(X)/dX is the same thing as d vec(A(X))/d vec(X).This is where the the relationship between the vec operator and Kronecker products is useful. Consider di?erentiating dx Ax with respect to A(rather than with respect to x as above):

d vec(x Ax) d vec(A)=

d(x ?x )vec(A)

d vec(A)

=(x ?x )

(the derivative of an m×n matrix A with respect to itself is I mn).

A more general product rule can be de?ned.Suppose that f: n→

m×p and g: n→ p×q,so f(x)g(x): n→ m×https://www.wendangku.net/doc/171169668.html,ing the relationship between the vec and Kronecker product operators

vec(I m f(x)g(x)I q)=(g(x) ?I m)vec(f(x))=(I q?f(x))vec(g(x)).

A natural product rule is therefore

Df(x)g(x)=(g(x) ?I m)f (x)+(I q?f(x))g (x).

This can be used to determine the derivative of dA A/dA where A is

m×n.

vec(A A)=(I n?A )vec(A)=(A ?I n)vec(A )=(A ?I n)T m,n vec(A).

Thus(using the product rule)

dA A

=(I n?A )+(A ?I n)T m,n.

This can be simpli?ed somewhat by noting that

(A ?I n)T m,n=T n,n(I n?A ).

Thus

dA A

=(I n2+T n,n)(I n?A ).

The product rule is also useful in determining the derivative of a matrix inverse:

dA?1A dA =(A ?I n)

dA?1

+(I n?A?1).

But A?1A is identically equal to I,so its derivative is identically0.Thus

dA?1

=?(A ?I n)?1(I n?A?1)=?(A? ?I n)(I n?A?1)=?(A? ?A?1).

It is also useful to have an expression for the derivative of a determinant. Suppose A is n×n with|A|=0.The determinant can be written as the

product of the i th row of the adjoint of A(A?)with the i th column of A: |A|=A?i·A·i.

Recall that the elements of the i th row of A?are not in?uenced by the elements in the i th column of A and hence

?|A|

?A·i

=A?i·.

To obtain the derivative with respect to all of the elements of A,we can concatenate the partial derivatives with respect to each column of A:

d|A| dA =[A?1·A?2·...A?n·]=|A|

[A?1]1·[A?1]2·...[A?1]n·

=|A|vec

The following result is an immediate consequence

d ln|A| dA =vec

Matrix di?erentiation results allow us to compute the derivatives of the solutions to certain classes of equilibrium problems.Consider,for example, the solution,x,to a linear complementarity problem LCP(M,q)that solves

Mx+q≥0,x≥0,x (Mx+q)=0

The i th element of x is either exactly equal to0or is equal to the i th element of Mx+q.De?ne a diagonal matrix D such that D ii=1if x>0and equal 0otherwise.The solution can then be written as x=??M?1Dq,where ?M=DM+I?D.It follows that

=??M?1D

and that

?x ?M =

??M?1

??M

=(?q D?I)(??M? ??M?1)(I?D) =q D?M? ??M?1D

=x ??x/?q

Given the prevalence of Kronecker products in matrix derivatives,it would be useful to have rules for computing derivatives of Kronecker prod-ucts themselves,i.e.dA?B/dA and dA?B/dB.Because each element of a Kronecker product involves the product of one element from A multiplied by one element of B,the derivative dA?B/dA must be composed of zeros and the elements of B arranged in a certain fashion.Similarly,the derivative dA?B/dB is composed of zeros and the elements of A arranged in a certain fashion.

It can be veri?ed that dA?B/dA can be written as

dA?B dA =

Ψ10 0

Ψ20 0

............

Ψq0 0

0Ψ1 0

0Ψ2 0

............

0Ψq 0

............

00...Ψ1

00...Ψ2

............

00...Ψq

=I n?

Ψ1

Ψ2

...

Ψq

where

Ψi=

B1i0 0

B2i0 0

............

B pi0 0

0B1i 0

0B2i 0

............

0B pi 0

............

00...B1i

00...B2i

............

00...B pi

=I m?B·i

This can be written more compactly as

dA?B

=(I n?T qm?I p)(I mn?vec(B))=(I nq?T mp)(I n?vec(B)?I m).

Similarly dA?B/dB can be written as

dA?B dB =

Θ10 0

0Θ1 0

............

00...Θ1

Θ20 0

0Θ2 0

............

00...Θ2

............

Θq0 0

0Θq 0

............

00...Θq

where

Θi=

A i10 0

0A i1 0

............

00...A i1

A i20 0

0A i2 0

............

00...A i2

............

A im0 0

0A im 0

............

00...A im

This can be written more compactly as

dA?B

=(I n?T qm?I p)(vec(A)?I pq)=(T pq?I mn)(I q?vec(A)?I p). Notice that if either A is a row vector(m=1)or B is a column vector

(q=1),the matrix(I n?T qm?I p)=Imnpq and hence can be ignored.

To illustrate a use for these relationships,consider the second derivative

of xx with respect to x,an n-vector.

dxx

=x?I n+I n?x;

hence

d2xx

dxdx

=(T nn?I n)(I n?vec(I n)+vec(I n)?I n).

Another example is

d2A?1 dAdA =(I n?T nn?I n)

I n2?vec

A?1

T nn+

vec

?I n2

A? ?A?1)

Often,especially in statistical applications,one encounters matrices that are symmetrical.It would not make sense to take a derivative with respect to the i,j th element of a symmetric matrix while holding the j,i th element constant.Generally it is preferable to work with a vectorized version of a symmetric matrix that excludes with the upper or lower portion of the matrix.The vech operator is typically taken to be the column-wise vector-ization with the upper portion excluded:

vech(A)=

A11

···

A n,1

A22

···

A n2

···

A nn?1

A nn

One obtains this by selecting elements of vec(A)and therefore we can write vech(A)=S n vec(A),

where S n is an n(n+1)/2×n2matrix of0s and1s,with a single1in each row.The vech operator can be applied to lower triangular matrices as well;there is no reason to take derivatives with respect to the upper part of a lower triangular matrix(it can also be applied to the transpose of an upper triangular matrix).The use of the vech operator is also important in e?cient computer storage of symmetric and triangular matrices.

To illustrate the use of the vech operator in matrix calculus applications, consider an n×n symmetric matrix C de?ned in terms of a lower triangular matrix,L,

C=LL .

Using the already familiar methods it is easy to see that

=(I n?L)T n,n+(L?I n).

Using the chain rule

d vech(C) d vech(L)=

d vech(C)

d vech(L)

=S n

S n.

Inverting this expression provides an expression for d vech(L)/d vech(C).

Related to matrix derivatives is the issue of Taylor expansions of matrix-to-matrix functions.One way to think of matrix derivatives is in terms of multidimensional arrays.An mn×pq matrix can also be thought of as an m×n×p×q4-dimensional array.The“reshape”function in MATLAB implements such transformations.The ordering of the individual elements has not change,only the way the elements are indexed.

The d th order Taylor expansion of a function f(X):R m×n→R p×q at ?X can be computed in the following way

f=vec(f(?X))

dX=vec(X??X)

for i=1to d

{

f i=f(i)(?X)dX

for j=2to i

f i=reshape(f i,mn(pq)j?2,pq)dX/j

f=f+f i

}

f=reshape(f,m,n)

The techniques described thus far can be applied to computation of derivatives of common“special”functions.First,consider the derivative of a nonnegative integer power A i of a square matrix A.Application of the chain rule leads to the recursive de?nition

dA i dA =

dA i?1A

=(A ?I)

dA i?1

+(I?A i?1)

which can also be expressed as a sum of i terms

dA i dA =

j=1

(A )i?j?A j?1.

This result can be used to derive an expression for the derivative of the matrix exponential function,which is de?ned in the usual way in terms of a Taylor expansion:

exp(A)=

∞

i=0

A i

Thus

d exp(A)

dA =

∞

i=0

j=1

(A )i?j?A j?1.

The same approach can be applied to the matrix natural logarithm:

ln(A)=?

∞

i=1

(I?A)i.

Summary of Operator Results

A is m×n,

B is p×q,X is de?ned comformably.

A?B=

A11B A12B...A1n B

A21B A22B...A2n B

............

A m1

B A m2B...A mn B

(AC?BD)=(A?B)(C?D) (A?B)?1=A?1?B?1

(A?B) =A ?B

vec(AXB)=(B ?A)vec(X) trace(AX)=vec(A ) vec(X) trace(AX)=trace(XA)

T m,n vec(A)=vec(A )

T n,m T m,n=I mn

T n,m=T?1

m,n

T m,n=T n,m

B?A=T p,m(A?B)T n,q

Summary of Di?erentiation Results

A is m×n,

B is p×q,x is n×1,X is de?ned comformably.

[Df]ij=d f i(x) dx j

dAx

D[αf(x)+βg(x)]=αDf(x)+βDg(x)

D[f(g(x))]=f (g(x))g (x)

D[f(x)g(x)]=(g(x) ?I m)f (x)+(I p?f(x))g (x). dx Ax

=x (A+A )

d vec(x Ax)

d vec(A)

=(x ?x )

dA A

=(I n2+T n,n)(I n?A )

dAA

=(I m2+T m,m)(A?I m)

dx A Ax

=2x ?x A (when x is a vector) dx AA x

=2x A?x (when x is a vector) dAXB

=B ?A

dA?1

=?(A? ?A?1)

d ln|A|

=vec(A? )

d trace(AX)

=vec(A )

dA?B

=(I n?T qm?I p)(I mn?vec(B))=(I nq?T mp)(I n?vec(B)?I m). dA?B

=(I n?T qm?I p)(vec(A)?I pq)=(T pq?I mn)(I q?vec(A)?I p).

dxx

=x?I n+I n?x

d2xx

dxdx

=(T nn?I n)(I n?vec(I n)+vec(I n)?I n)

dA i dA =(A ?I)

dA i?1

+(I?A i?1)=

j=1

(A )i?j?A j?1.

d exp(A)

dA =

∞

i=0

j=1

(A )i?j?A j?1.

实现稀疏矩阵(采用三元组表示)的基本运算实验报告

实现稀疏矩阵（采用三元组表示）的基本运算实验报告一实验题目: 实现稀疏矩阵（采用三元组表示）的基本运算二实验要求: （1）生成如下两个稀疏矩阵的三元组 a 和 b；（上机实验指导 P92 ）（2）输出 a 转置矩阵的三元组；（3）输出a + b 的三元组；（4）输出 a * b 的三元组；三实验内容: 稀疏矩阵的抽象数据类型: ADT SparseMatrix { 数据对象:D={aij| i = 1,2,3,….,m; j =1,2,3,……,n; ai,j∈ElemSet,m和n分别称为矩阵的行数和列数} 数据关系: R={ Row , Col } Row ={ | 1≤i≤m , 1≤j ≤n-1} Col ={| 1≤i≤m-1,1≤j ≤n} 基本操作:

CreateSMatrix(&M) 操作结果：创建稀疏矩阵 M PrintSMatrix(M) 初始条件：稀疏矩阵M已经存在操作结果：打印矩阵M DestroySMatrix(&M) 初始条件：稀疏矩阵M已经存在操作结果：销毁矩阵M CopySMatrix(M, &T) 初始条件：稀疏矩阵M已经存在操作结果：复制矩阵M到T AddSMatrix(M, N, &Q) 初始条件：稀疏矩阵M、N已经存在操作结果：求矩阵的和Q=M+N SubSMatrix(M, N, &Q) 初始条件：稀疏矩阵M、N已经存在操作结果：求矩阵的差Q=M-N TransposeSMatrix(M, & T) 初始条件：稀疏矩阵M已经存在操作结果：求矩阵M的转置T MultSMatrix(M, N, &Q) 初始条件：稀疏矩阵M已经存在