SIERE a hybrid semi-implicit exponential integrator for efficiently simulating stiff deformable objects
3 METHOD
3.1 Model reduction and subspace splitting
Next, we define the splitting $G$ and $H$, crucial to the success of our method. The idea is to apply ERE in the subspace of the first s modes ($s ≪ n$: typically, $5 ≤ s ≤ 20$) and project it back to the original full space. In the bridge example of Figure 1, $n ≈ 100, 000$, so this is a rather large reduction.
Then we use SI on the remaining unevaluated part, as per Eq. (5). We use mass-PCA to find our reduced space. That is, considering at the beginning of each time step a solution mode of the form $q(t) = w exp(ı\sqrt{λ}t)$ for the ODE ${\proselabel{siere}{{M}}} q + {\proselabel{siere}{{K}}} q = 0$, we solve the generalized eigenvalue problem
$$ {\proselabel{siere}{{K}}} w = λ {\proselabel{siere}{{M}}} w \notag$$
for the $s$ smallest eigenvalues $λ$ and their corresponding eigenvectors $w$ (dominant modes). In Matlab, for instance, this can be achieved by calling the function eigs. In our implementation, we use the C++ Spectra library [Qiu 2019]. Denote this partial spectral decomposition by
$${\proselabel{siere}{{K}}} {\proselabel{siere}{{U_s}}} ={\proselabel{siere}{{M}}} {\proselabel{siere}{{U_s}}} Λ_s\tag{6}\label{6}$$
where the “long and skinny” ${\proselabel{siere}{{U_s}}}$ is $n × s$, has the first s eigenvectors $w$ as its columns, and the small $Λ_s$ is a diagonal $s × s$ matrix with the eigenvalues $λ_1, ..., λ_s$ on the diagonal. Notice that both ${\proselabel{siere}{{K}}}$ and ${\proselabel{siere}{{M}}}$ are large sparse symmetric matrices. In addition, ${\proselabel{siere}{{M}}}$ is positive definite, so ${\proselabel{siere}{{U_s}}}$ has M-orthogonal columns:
$${\proselabel{siere}{{U_s}}} ^T {\proselabel{siere}{{M}}} {\proselabel{siere}{{U_s}}} = I_s, {\proselabel{siere}{{U_s}}} ^T {\proselabel{siere}{{K}}} {\proselabel{siere}{{U_s}}} = Λ_s. \tag{7}\label{7}$$
Next, we write Eq. (2) in the split form Eq. (3), with the splitting $H$ and $G$ defined based on the partial spectral decomposition Eq. (6). We define at each time step
$$\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\begin{align*}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$G(u)$', 'siere', 'def', false, '')", "id":"siere-$G(u)$", "sym":"$G(u)$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {G(u)} } & = \begin{bmatrix}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$v_G$', 'siere', 'use', false, '')", "id":"siere-$v_G$", "sym":"$v_G$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {v_G} }\\
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'M', 'siere', 'use', false, '')", "id":"siere-M", "sym":"M", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{M}} }^{-1}\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$f_G$', 'siere', 'use', false, '')", "id":"siere-$f_G$", "sym":"$f_G$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {f_G} }\\
\end{bmatrix}\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['$f_G$', '$v_G$', 'M', '$G(u)$'], false, [], [], 'YCRHKHUpJGAgPSBbYCR2X0ckYAogICAgICAgICAgTeKBu8K5YCRmX0ckYF0=');"} }{}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$H(u)$', 'siere', 'def', false, '')", "id":"siere-$H(u)$", "sym":"$H(u)$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {H(u)} } & = \begin{bmatrix}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$v_H$', 'siere', 'use', false, '')", "id":"siere-$v_H$", "sym":"$v_H$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {v_H} }\\
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'M', 'siere', 'use', false, '')", "id":"siere-M", "sym":"M", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{M}} }^{-1}\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$f_H$', 'siere', 'use', false, '')", "id":"siere-$f_H$", "sym":"$f_H$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {f_H} }\\
\end{bmatrix}\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['$f_H$', '$v_H$', 'M', '$H(u)$'], false, [], [], 'YCRIKHUpJGAgPSBbYCR2X0gkYAogICAgICAgICAgTeKBu8K5YCRmX0gkYF0=');"} }{}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$v_G$', 'siere', 'def', false, '')", "id":"siere-$v_G$", "sym":"$v_G$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {v_G} } & = \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }}^T\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'M', 'siere', 'use', false, '')", "id":"siere-M", "sym":"M", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{M}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'v', 'siere', 'use', false, '')", "id":"siere-v", "sym":"v", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{v}} }\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['v', '$U_s$', 'M', '$v_G$'], false, [], [], 'YCR2X0ckYCA9IGAkVV9zJGBgJFVfcyRgXlQgTXY=');"} }{}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$v_H$', 'siere', 'def', false, '')", "id":"siere-$v_H$", "sym":"$v_H$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {v_H} } & = \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'v', 'siere', 'use', false, '')", "id":"siere-v", "sym":"v", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{v}} } - \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$v_G$', 'siere', 'use', false, '')", "id":"siere-$v_G$", "sym":"$v_G$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {v_G} }\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['$v_G$', 'v', '$v_H$'], false, [], [], 'YCR2X0gkYCA9ICB2IC0gYCR2X0ckYA==');"} }{}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$f_G$', 'siere', 'def', false, '')", "id":"siere-$f_G$", "sym":"$f_G$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {f_G} } & = \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'M', 'siere', 'use', false, '')", "id":"siere-M", "sym":"M", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{M}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }}^T\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'f', 'siere', 'use', false, '')", "id":"siere-f", "sym":"f", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{f}} }\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['f', '$U_s$', 'M', '$f_G$'], false, [], [], 'YCRmX0ckYCA9ICBNYCRVX3MkYGAkVV9zJGBeVCBm');"} }{}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$f_H$', 'siere', 'def', false, '')", "id":"siere-$f_H$", "sym":"$f_H$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {f_H} } & = \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'f', 'siere', 'use', false, '')", "id":"siere-f", "sym":"f", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{f}} } - \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$f_G$', 'siere', 'use', false, '')", "id":"siere-$f_G$", "sym":"$f_G$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {f_G} }\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['$f_G$', 'f', '$f_H$'], false, [], [], 'YCRmX0gkYCA9ICBmIC0gYCRmX0ckYA==');"} }{}
\end{align*}
\tag{8}\label{8}$$
We also need the Jacobian matrices
$$\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\begin{align*}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$J_G$', 'siere', 'def', false, '')", "id":"siere-$J_G$", "sym":"$J_G$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {J_G} } & = \begin{bmatrix}
0 & \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }}^T\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'M', 'siere', 'use', false, '')", "id":"siere-M", "sym":"M", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{M}} }\\
-\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }}^T\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'K', 'siere', 'use', false, '')", "id":"siere-K", "sym":"K", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{K}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }}^T\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'M', 'siere', 'use', false, '')", "id":"siere-M", "sym":"M", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{M}} } & 0\\
\end{bmatrix}\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['$U_s$', 'M', 'K', '$J_G$'], false, [], [], 'YCRKX0ckYCA9IFswICAgIGAkVV9zJGBgJFVfcyRgXlRNCiAgICAgIC1gJFVfcyRgYCRVX3MkYF5US2AkVV9zJGBgJFVfcyRgXlRNIDAgXQ==');"} }{}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$J_H$', 'siere', 'def', false, '')", "id":"siere-$J_H$", "sym":"$J_H$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {J_H} } & = \begin{bmatrix}
0 & I_{ \mathit{n} }\\
-\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'M', 'siere', 'use', false, '')", "id":"siere-M", "sym":"M", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{M}} }^{-1}\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'K', 'siere', 'use', false, '')", "id":"siere-K", "sym":"K", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{K}} } & 0\\
\end{bmatrix} - \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$J_G$', 'siere', 'use', false, '')", "id":"siere-$J_G$", "sym":"$J_G$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {J_G} }\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['$J_G$', 'M', 'K', '$J_H$'], false, [], [], 'YCRKX0gkYCA9ICBbMCAgICAgSV9uCiAgICAgICAgICAgIC1N4oG7wrlLIDBdIC0gYCRKX0ckYA==');"} }{}
\end{align*}
\tag{9}\label{9}$$
Notice that the ERE expression, ${\proselabel{siere}{{h}}} {\proselabel{siere}{{φ_1}}} ( {\proselabel{siere}{{h}}} {\proselabel{siere}{{J_G}}} ) {\proselabel{siere}{{G(u)}}} $, can be evaluated in the subspace first, and then projected back to the original space.
The additive method defined by inserting Eqs. (8) and (9) into Eq.(5) has three advantages:
- (1) At each time step, the majority of the update comes from ERE in the dominating modes. Thus it is less affected by artificial damping from SI.
- (2) The computation load of ERE is greatly reduced, because the stiff part is handled by SI (or BE for that matter). Furthermore, the evaluation of the exponential function in the subspace has only marginal cost since the crucial matrix involved has been diagonalized.
- (3) The “warm start” for SI makes its result closer to that of BE.
ERE update in the subspace: To evaluate the update in the subspace of dimension s we rewrite Eq. (5) as
$$\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\begin{align*}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$u_+$', 'siere', 'def', false, '')", "id":"siere-$u_+$", "sym":"$u_+$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {u_+} } & = \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'u', 'siere', 'use', false, '')", "id":"siere-u", "sym":"u", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{u}} } + \left( \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$\\\\boldsymbol{I}$', 'siere', 'use', false, '')", "id":"siere-$\\\\boldsymbol{I}$", "sym":"$\\\\boldsymbol{I}$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\boldsymbol{I}} } - \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'h', 'siere', 'use', false, '')", "id":"siere-h", "sym":"h", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{h}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$J_H$', 'siere', 'use', false, '')", "id":"siere-$J_H$", "sym":"$J_H$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {J_H} } \right)^{-1}\left( \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'h', 'siere', 'use', false, '')", "id":"siere-h", "sym":"h", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{h}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$H(u)$', 'siere', 'use', false, '')", "id":"siere-$H(u)$", "sym":"$H(u)$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {H(u)} } + \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'h', 'siere', 'use', false, '')", "id":"siere-h", "sym":"h", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{h}} }\begin{bmatrix}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} } & 0\\
0 & \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }\\
\end{bmatrix}\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$φ_1$', 'siere', 'use', false, '')", "id":"siere-$φ_1$", "sym":"$φ_1$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {φ_1} }\left( \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'h', 'siere', 'use', false, '')", "id":"siere-h", "sym":"h", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{h}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$J_G^r$', 'siere', 'use', false, '')", "id":"siere-$J_G^r$", "sym":"$J_G^r$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {J_G^r} } \right)\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$G^r(u)$', 'siere', 'use', false, '')", "id":"siere-$G^r(u)$", "sym":"$G^r(u)$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {G^r(u)} } \right)\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['$J_G^r$', '$φ_1$', '$\\\\boldsymbol{I}$', '$J_H$', '$U_s$', '$H(u)$', '$G^r(u)$', 'h', 'u', '$u_+$'], false, [], [], 'YCR1XyskYCA9ICB1ICsgKGAkXGJvbGRzeW1ib2x7SX0kYCAtaGAkSl9IJGAp4oG7wrkoaCBgJEgodSkkYCArIGhbYCRVX3MkYCAwCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgMCAgIGAkVV9zJGBdIGAkz4ZfMSRgKGhgJEpfR15yJGApIGAkR15yKHUpJGAp');"} }{}
\end{align*}
\tag{10}\label{10}$$
where
$$\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\begin{align*}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$J_G^r$', 'siere', 'def', false, '')", "id":"siere-$J_G^r$", "sym":"$J_G^r$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {J_G^r} } & = \begin{bmatrix}
0 & I_{ \mathit{s} }\\
-{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }}^T\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'K', 'siere', 'use', false, '')", "id":"siere-K", "sym":"K", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{K}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} } & 0\\
\end{bmatrix}\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['$U_s$', 'K', '$J_G^r$'], false, [], [], 'YCRKX0deciRgID0gWzAgICAgIElfcwogICAgICAgICAgICAtYCRVX3MkYF5US2AkVV9zJGAgMF0=');"} }{}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$G^r(u)$', 'siere', 'def', false, '')", "id":"siere-$G^r(u)$", "sym":"$G^r(u)$", "func":"siere", "localFunc":"", "type":"def", "case":"equation"} }{ {G^r(u)} } & = \begin{bmatrix}
{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }}^T\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'M', 'siere', 'use', false, '')", "id":"siere-M", "sym":"M", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{M}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'v', 'siere', 'use', false, '')", "id":"siere-v", "sym":"v", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{v}} }\\
{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'siere', 'use', false, '')", "id":"siere-$U_s$", "sym":"$U_s$", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }}^T\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'f', 'siere', 'use', false, '')", "id":"siere-f", "sym":"f", "func":"siere", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{f}} }\\
\end{bmatrix}\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'siere', ['v', 'f', '$U_s$', 'M', '$G^r(u)$'], false, [], [], 'YCRHXnIodSkkYCA9IFtgJFVfcyRgXlRNdgogICAgICAgICAgICBgJFVfcyRgXlRmXQ==');"} }{}
\end{align*}
\tag{11}\label{11}$$
The evaluation of the action of the matrix function ${\proselabel{siere}{{φ_1}}}$ involves only matrices of size $2s × 2s$. Furthermore, the matrix function ${\proselabel{siere}{{φ_1}}}$ can be evaluated directly through the eigenpairs of ${\proselabel{siere}{{J_G^r}}}$
$$
\Bigg\{ ı\sqrt{λ_l}, \begin{bmatrix} e_l \\ ı\sqrt{λ_l}e_l \end{bmatrix} \Bigg\} ,
\Bigg\{ -ı\sqrt{λ_l}, \begin{bmatrix} -e_l \\ ı\sqrt{λ_l}e_l \end{bmatrix} \Bigg\},
l = 1,...,s,
\tag{12}\label{12}$$
with $e_l$ being the $l^{th}$ column of the identity matrix.
The large $n×n$ linear system solved in Eq. (10) is not sparse due to the fill-in introduced by the small subspace projection. Specifically, the off-diagonal blocks of the Jacobian matrix ${\proselabel{siere}{{J_G}}}$ defined in Eq. (9a) are not sparse. If not treated carefully, solving the linear system in Eq. (5) and Eq. (10) can be prohibitively costly. Fortunately, this modification matrix has the low rank $s$. We can write
$$\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\begin{align*}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$J_G$', 'second', 'def', false, '')", "id":"second-$J_G$", "sym":"$J_G$", "func":"second", "localFunc":"", "type":"def", "case":"equation"} }{ {J_G} } & = \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$Y_1$', 'second', 'use', false, '')", "id":"second-$Y_1$", "sym":"$Y_1$", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {Y_1} }{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$Z_1$', 'second', 'use', false, '')", "id":"second-$Z_1$", "sym":"$Z_1$", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {Z_1} }}^T + \idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$Y_2$', 'second', 'use', false, '')", "id":"second-$Y_2$", "sym":"$Y_2$", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {Y_2} }{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$Z_2$', 'second', 'use', false, '')", "id":"second-$Z_2$", "sym":"$Z_2$", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {Z_2} }}^T\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'second', ['$Y_2$', '$Y_1$', '$Z_1$', '$Z_2$', '$J_G$'], false, [], [], 'YCRKX0ckYCA9IGAkWV8xJGBgJFpfMSRgXlQgKyBgJFlfMiRgYCRaXzIkYF5U');"} }{}
\end{align*}
\notag$$
where
$$\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\begin{align*}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$Y_1$', 'second', 'def', false, '')", "id":"second-$Y_1$", "sym":"$Y_1$", "func":"second", "localFunc":"", "type":"def", "case":"equation"} }{ {Y_1} } & = \begin{bmatrix}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'second', 'use', false, '')", "id":"second-$U_s$", "sym":"$U_s$", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }\\
0\\
\end{bmatrix}\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'second', ['$U_s$', '$Y_1$'], false, [], [], 'YCRZXzEkYCA9ICBbYCRVX3MkYAogICAgICAgICAgICAgMF0=');"} }{}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$Z_1$', 'second', 'def', false, '')", "id":"second-$Z_1$", "sym":"$Z_1$", "func":"second", "localFunc":"", "type":"def", "case":"equation"} }{ {Z_1} } & = \begin{bmatrix}
0\\
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'M', 'second', 'use', false, '')", "id":"second-M", "sym":"M", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{M}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'second', 'use', false, '')", "id":"second-$U_s$", "sym":"$U_s$", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }\\
\end{bmatrix}\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'second', ['$U_s$', 'M', '$Z_1$'], false, [], [], 'YCRaXzEkYCA9ICBbIDAKICAgICAgICAgICAgIE1gJFVfcyRgXQ==');"} }{}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$Y_2$', 'second', 'def', false, '')", "id":"second-$Y_2$", "sym":"$Y_2$", "func":"second", "localFunc":"", "type":"def", "case":"equation"} }{ {Y_2} } & = \begin{bmatrix}
0\\
-\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'second', 'use', false, '')", "id":"second-$U_s$", "sym":"$U_s$", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }{\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'second', 'use', false, '')", "id":"second-$U_s$", "sym":"$U_s$", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }}^T\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'K', 'second', 'use', false, '')", "id":"second-K", "sym":"K", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{K}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'second', 'use', false, '')", "id":"second-$U_s$", "sym":"$U_s$", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }\\
\end{bmatrix}\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'second', ['$U_s$', 'K', '$Y_2$'], false, [], [], 'YCRZXzIkYCA9ICBbMAogICAgICAgICAgICAtYCRVX3MkYGAkVV9zJGBeVEtgJFVfcyRgXQ==');"} }{}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$Z_2$', 'second', 'def', false, '')", "id":"second-$Z_2$", "sym":"$Z_2$", "func":"second", "localFunc":"", "type":"def", "case":"equation"} }{ {Z_2} } & = \begin{bmatrix}
\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, 'M', 'second', 'use', false, '')", "id":"second-M", "sym":"M", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {\mathit{M}} }\idlabel{ {"onclick":"event.stopPropagation(); onClickSymbol(this, '$U_s$', 'second', 'use', false, '')", "id":"second-$U_s$", "sym":"$U_s$", "func":"second", "localFunc":"", "type":"use", "case":"equation"} }{ {U_s} }\\
0\\
\end{bmatrix}\\\eqlabel{ {"onclick":"event.stopPropagation(); onClickEq(this, 'second', ['$U_s$', 'M', '$Z_2$'], false, [], [], 'YCRaXzIkYCA9ICBbIE1gJFVfcyRgCiAgICAgICAgICAgICAwXQ==');"} }{}
\end{align*}
\notag$$
The linear system in Eq. (10) becomes
$$ I−h J_H =(I−h J)+h {\proselabel{second}{{Y_1}}} {\proselabel{second}{{Z_1}}} ^T +h {\proselabel{second}{{Y_2}}} {\proselabel{second}{{Z_2}}} ^T \tag{13}\label{13}$$
where the four matrices $Y_i$ and $Z_i$ are all “long and skinny” like ${\proselabel{second}{{U_s}}}$, while the matrix $J$ is square and large, but very sparse. Figure 3 illustrates this situation. For the linear system to be solved in Eq. (10) we may employ an iterative method such as conjugate gradient, whereby the matrix-vector products involving $J$ or $Y_iZ_i^T$ are all straightforward to carry out efficiently. However, we have often found out that a direct solution method is more appropriate for these linear equations in our context. In our implementation we use pardiso [De Coninck et al. 2016; Kourounis et al. 2018; Verbosio et al. 2017]. For this we can employ the formula of Sherman, Morrison and Woodbury (SMW) [Nocedal and Wright 2006], given by
$$ (A+YZ^T)^{−1} =A^{−1}−A^{−1}Y(I +Z^TA^{−1}Y)^{−1}Z^TA^{−1} \tag{14}\label{14}$$
to solve the linear system in Eq. (10). In our specific notation we set at each time step $A = I − h J_H$ in Eq. (14), and apply the formula twice:once for $Y = {\proselabel{second}{{Y_1}}}$, $Z = {\proselabel{second}{{Z_1}}}$, and once for $Y = {\proselabel{second}{{Y_2}}}$, $Z = {\proselabel{second}{{Z_2}}}$. Note that the matrices $I + Z^T A^{−1}Y$ in Eq. (14) are only $2s × 2s$, and this results in an efficient implementation, so long as the subspace dimension s remains small.