суммировать матричный ряд: как увеличить цикл матричного умножения и экспоненциального

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: Уценка (требуется поддержка MathJax) для рисунка 1

##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: Уценка (нужна поддержка MathJax) для картинки 2

$\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}.$$