对矩阵系列求和：如何提升矩阵乘法和指数循环

Question

我有一个优化问题，我使用optim 例程和卡尔曼滤波器 来找到最大似然参数。这是一个处理矩阵的问题，由于随机性，我必须使用积分来估计基础状态变量的矩阵。

我在 R 中找不到处理矩阵积分的合适函数。如果我自己编程它会相对较慢，从而使优化算法更慢。我最感兴趣的是速度。

我尝试了以下方法。我生成了随机矩阵Kappa和sigma_mat，expm是一个计算矩阵指数的函数。如何使这个循环更快？当然我减少了迭代次数，但我想在估计积分时保持一些准确性。

install.packages("expm")
library(expm) #This loads the function to calculate matrix exponentials

#Generate some example random matrices
set.seed(0)
Kappa <- matrix(rnorm(9), nrow = 3, ncol = 3)
sigma_mat <- matrix(rnorm(9), nrow = 3, ncol = 3)

#Now we estimate the integral
Q <- 0
for(i in seq(0,1,length.out=5000)*(1/12)){
  Q <- Q + expm(-Kappa*i) %*% sigma_mat %*% t(sigma_mat) %*% expm(-t(Kappa)*i)
  }
Q <- Q / 5000
Q

Answer 1

您的原始代码和速度

set.seed(0)
library(expm)
Kappa <- matrix(rnorm(9), 3, 3)
sigma_mat <- matrix(rnorm(9), 3, 3)

system.time({
Q1 <- 0
for(i in seq(0,1,length.out=5000)*(1/12)){
  Q1 <- Q1 + expm(-Kappa*i) %*% sigma_mat %*% t(sigma_mat) %*% expm(-t(Kappa)*i)
  }
Q1 <- Q1 / 5000
Q1})
#   user  system elapsed 
#  4.464   0.136   4.605 

---

---

以下 R 代码实现了上述算法。变量名与上述推导中使用的一致。

system.time({
A <- -Kappa
B <- sigma_mat
E <- eigen(expm(A))
d <- E[[1]]
U <- E[[2]]
C <- tcrossprod(solve(U, B))
K <- tcrossprod(d)
a <- 0
b <- 1 / 12
n <- 5000
W <- K ^ {(b - a) / (n - 1)}
Q2 <- (1 - W ^ n) / (1 - W)
Q2 <- C * Q2
Q2 <- Re(tcrossprod(U %*% Q2, U))
Q2 <- Q2 / n
Q2})
#   user  system elapsed 
#  0.004   0.000   0.002 

## check that the computational result is correct
all.equal(Q1, Q2)
#[1] TRUE

---

附录 1：图片 1 的 Markdown（需要 MathJax 支持）

##Notation

Let $X$ be a square matrix,

 - $X'$ is the transpose of $X$;
 - $X^{-1}$ is the inverse of $X$;
 - $X^i$ is the i-th power of $X$. For example, $X^3 = X * X * X$ where $*$ is the matrix multiplication;
 - $X^{[i]}$ is the element-wise i-th power of $X$. For example, $X^{[3]} = X \circ X \circ X$ where $\circ$ is Hadamard matrix product, i.e., element-wise matrix product;
 - $\exp(X)$ is the matrix exponential;
 - $\texttt{diag}(X)$ is the main diagonal vector of $X$.

Note that both matrix power and its element-wise version can be defined for non-integer $i$ value.

----

##Mathematical derivation

Let matrix $A$ be your `-Kappa` and matrix $B$ be your `sigma_mat`, you are computing $$\sum_i\exp(Ai)BB'\exp(A'i) = \sum_i\big(\exp(Ai)B\big)\big(\exp(Ai)B\big)' = \sum_i\big(\exp(A)^iB\big)\big(\exp(A)^iB\big)'.$$ Consider an eigen decomposition $\exp(A) = UDU^{-1}$, then the summation becomes $$\sum_i\big(UD^iU^{-1}B\big)\big(UD^iU^{-1}B\big)' = U\big(\sum_iD^iCD^i\big)U' = U\big(C \circ \sum_i K^{[i]}\big)U',$$ where $C = \big(U^{-1}B\big)\big(U^{-1}B\big)'$, $K = dd'$ and $d = \texttt{diag}(D)$.

附录 2：图片 2 的 Markdown（需要 MathJax 支持）

$\sum_iK^{[i]}$ is the sum of a geometric series and an analytical solution exists. Suppose $i$ takes $n$ evenly spaced values on $\left[a, b\right]$, that is, $i = a,\ a + j,\ a + 2j,\ \cdots,\ a + (n - 1)j,$ where $j = \frac{b - a}{(n-1)}$. Let $W = K^{[j]}$, there is $$\sum_iK^{[i]} = K^{[a]}\circ\sum_{j = 0}^{n - 1}W^{[j]} = K^{[a]} \circ \frac{1 - W^{[n]}}{1 - W}.$$