Derivative of a Scalar Function with respect to a Matrix

Description

This is the wrapper function that calls mD0 function. It pre-processes the input and post-processes the output .

Usage

mD(
  expr,
  X_ = "X",
  sym = 0,
  simplify = 1,
  trace_chain = 1,
  print = 0,
  debug = 0
)

Arguments

Details

To numerically check the result, use
gramdmn and gradma
functions. (See example below.)

Currently, tr and det are the only functions available as the input.

Whe simplify=1 is given,
Expand, distribute_t and remove_paren
will be used to simplify the result.

Value

a string containing the first derivative of rss w.r.t. X_.

Examples

# we need this.
modify_math_operators()

mD("tr(X)","X")
mD("tr(A%*%X)","X")
mD("tr(A%*%t(X))","X")
mD("tr(X%*%A)","X")
mD("tr(t(X)%*%A)","X")
mD("tr(A%*%t(X)%*%B%*%X%*%C)","X")
mD("tr(A%*%X%*%B%*%X%*%C)","X")

mD("tr(A%*%inv(X))","X")
mD("tr(A%*%inv(t(X)))","X")
mD("tr(A%*%inv(B%*%X%*%C)%*%D)","X")
mD("tr(A%*%inv(t(X)%*%B%*%X)%*%C)","X")

mD("tr(Diag(t(X)%*%X))","X")
mD("tr(A%*%(Diag(t(X)%*%C%*%X))%*%B)","X")

mD("tr((X%.%B))","X")
mD("tr(A%*%(X%.%B)%*%C)","X")
mD("tr((t(X)%.%X))","X")
mD("tr(A%*%(t(X)%.%B%.%X)%*%C)","X")

mD("det(X)","X")
mD("det(t(X))","X")
mD("det(A%*%X)","X")

# symmetric X
mD("tr(A%*%X)", "X", sym=1)
mD("tr(t(X)%*%A%*%X)", "X", sym=1)



## Not run: 
 # Must copy and paste to run the following examples.

 # Examples of numerical checking

 # inverse: tr(A%*%inv(t(X)%*%B%*%X)%*%C)
 expr <- "tr(A%*%inv(t(X)%*%B%*%X)%*%C)"
 dexpr <- mD(expr,"X")
 printm(dexpr)
 nvar1 <- 2; nvar2 <- 3
 set.seed(1701)
 Xn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 An <- matrix(rnorm(nvar2*nvar1),nvar2,nvar1)
 Bn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 Cn <- matrix(rnorm(nvar1*nvar2),nvar1,nvar2)
 values <- "X=Xn; A=An; B=Bn; C=Cn"
 Ga <- 0; Gn <- 9999
 Ga <- gradma( expr, dexpr=dexpr, values=values)
 Gn <- gradmn( expr, values=values)
 printm(Ga,Gn)
 printm(max(abs(Ga-Gn))/mean(abs(Gn)))


 # diag: tr(A%*%(Diag(t(X)%*%C%*%X))%*%B)
 expr <- "tr(A%*%(Diag(t(X)%*%C%*%X))%*%B)"
 dexpr <- mD(expr, "X")
 printm(dexpr)
 nvar1 <- 2; nvar2 <- 3
 set.seed(1701)
 Xn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 An <- matrix(rnorm(nvar2*nvar1),nvar2,nvar1)
 Bn <- matrix(rnorm(nvar1*nvar2),nvar1,nvar2)
 Cn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 values <- "X=Xn; A=An; B=Bn; C=Cn; I=In"
 rm(Ga,Gn)
 Ga=Eval(dexpr, values=values, fullsymb=1)
 Gn=gradmn( expr, values=values )
 printm(Ga,Gn)
 printm(max(abs(Ga-Gn))/mean(abs(Gn)))


 # Hadamar product: tr(A%*%(t(X)%.%C%.%X)%*%B)
 expr <- "tr(A%*%(t(X)%.%C%.%X)%*%B)"
 dexpr <- mD(expr, "X")
 printm(dexpr)
 nvar1 <- 2; nvar2 <- 3
 set.seed(1701)
 Xn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 An <- matrix(rnorm(nvar2*nvar1),nvar2,nvar1)
 Bn <- matrix(rnorm(nvar1*nvar2),nvar1,nvar2)
 Cn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 values <- "X=Xn; A=An; B=Bn; C=Cn; D=Dn"
 rm(Ga,Gn)
 Ga=Eval(dexpr, values=values, fullsymb=1)
 Gn=gradmn( expr, values=values )
 printm(Ga,Gn)
 printm(max(abs(Ga-Gn))/mean(abs(Gn)))


 # determinant: det(A%*%X%*%B)+det(C%*%t(X)%*%D)
 expr <- "det(A%*%X%*%B)+det(C%*%t(X)%*%D)"
 dexpr <- mD(expr,"X")
 printm(dexpr)
 nvar1 <- 3
 set.seed(1701)
 Xn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 An <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 Bn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 Cn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 Dn <- matrix(rnorm(nvar1*nvar1),nvar1,nvar1)
 values <- "X=Xn; A=An; B=Bn; C=Cn; D=Dn"
 Ga <- 0; Gn <- 9999
 Ga <- gradma( expr, dexpr=dexpr, values=values)
 Gn <- gradmn( expr, values=values)
 printm(Ga,Gn)
 printm(max(abs(Ga-Gn))/mean(abs(Gn)))



# Some complex examples

# rss of doubly weighted multivariate multiple regression
rss <- ssq("Y-X%*%B", U="W", V="V", expand=2) |> printm()

# differantiate rss w.r.t the regression coefficient matrix B.
dexpr <- mD(rss, "B", trace_chain=1, print=1)
printm(dexpr)
# numerical checking by gradma and gradmn functions
set.seed(1701)
n <- 5; ny <- 4; nx <- 3
Vn <- diag(ny)*2
Wn <- diag(n)*3
Xn <- matrix(rnorm(n*nx),n,nx)
Bn <- matrix(rnorm(nx*ny),nx,ny)
Yn <- matrix(rnorm(n*ny),n,ny)
values <- c("B=Bn; Y=Yn; X=Xn; V=Vn; W=Wn")#'
Ga <- 0; Gn <- 9999
Ga <- gradma( rss, dexpr=dexpr, values=values)
Gn <- gradmn( rss, values=values)
printm(Ga,Gn)
printm(max(abs(Ga-Gn))/mean(abs(Gn)))


# LS fa
rss <- ssq("S-(L%*%Phi%*%t(L)+Psi)", expand=2)

# d ssq / d L
dexpr <- mD(rss,"L", trace_chain=1, print=1)
printm(dexpr)
set.seed(1701)
nvar <- 3; ndim <- 2
Ln <- matrix(rnorm(nvar*ndim),nvar,ndim)
Phin <- matrix(runif(ndim*ndim),ndim,ndim)
Phin <- (Phin+t(Phin))/2; Phin <- Phin-Diag(Phin)+diag(ndim)
Psin <- diag(nvar) / 10
Sn <- matrix(runif(nvar*nvar),nvar,nvar); Sn=(Sn+t(Sn))/2; diag(Sn)=1
values <- "L=Ln; Phi=Phin; Psi=Psin; S=Sn"
Ga <- 0; Gn <- 9999
Ga <- gradma( rss, dexpr=dexpr, values=values)
Gn <- gradmn( rss, values=values)
printm(Ga,Gn)
printm(max(abs(Ga-Gn))/mean(abs(Gn)))

# d ssq / d Phi
dexpr <- mD(rss,"Phi", print=1)
printm(dexpr)
values <- "Phi=Phin; L=Ln; Psi=Psin; S=Sn"
Ga <- 0; Gn=9999
Ga <- gradma( rss, dexpr=dexpr, values=values)
Gn <- gradmn( rss, values=values)
printm(Ga,Gn)
printm(max(abs(Ga-Gn))/mean(abs(Gn)))

# d ssq / d Psi
dexpr <- mD(rss,"Psi", print=1)
printm(dexpr)
values <- "Psi=Psin; Phi=Phin; L=Ln; S=Sn"
Ga <- 0; Gn <- 9999
Ga <- gradma( rss, dexpr=dexpr, values=values)
Gn <- gradmn( rss, values=values)
printm(Ga,Gn)
printm(max(abs(Ga-Gn))/mean(abs(Gn)))



# ML fa
expr1 <- "tr(inv(L%*%Phi%*%t(L)+Psi)%*%S)"
expr2 <- "log(det(inv(L%*%Phi%*%t(L)+Psi)%*%S))"
F <- expr1+expr2

# d F / d L
dexpr <- mD(F,"L")
printm(dexpr)
set.seed(1701)
nvar <- 3; ndim <- 2
Ln <- matrix(rnorm(nvar*ndim),nvar,ndim)
Phin <- matrix(runif(ndim*ndim),ndim,ndim)
Phin <- (Phin+t(Phin))/2; Phin <- Phin-Diag(Phin)+diag(ndim)
Psin <- diag(nvar) / 10
Sn <- matrix(runif(nvar*nvar),nvar,nvar); Sn=(Sn+t(Sn))/2; diag(Sn)=1
values <- "L=Ln; Phi=Phin; Psi=Psin; S=Sn"
Ga <- 0; Gn <- 9999
Ga <- gradma( F, dexpr=dexpr, values=values)
Gn <- gradmn( F, values=values)
printm(Ga,Gn)
printm(max(abs(Ga-Gn))/mean(abs(Gn)))


# d F / d Phi
dexpr <- mD(F,"Phi", sym=1)
printm(dexpr)
values <- "Phi=Phin; L=Ln; Psi=Psin; S=Sn"
Ga <- 0; G <- 9999
Ga <- gradma( F, dexpr=dexpr, values=values)
Gn <- gradmn( F, values=values, sym=1)
printm(Ga,Gn)
printm(max(abs(Ga-Gn))/mean(abs(Gn)))


# d F / d Psi
dexpr <- mD(F,"Psi", sym=1)
printm(dexpr)
values <- "Psi=Psin; Phi=Phin; L=Ln; S=Sn"
Ga <- 0; G <- 9999
Ga <- gradma( F, dexpr=dexpr, values=values)
Gn <- gradmn( F, values=values, sym=1)
printm(Ga,Gn)
printm(max(abs(Ga-Gn))/mean(abs(Gn)))


## End(Not run)