When solving a system of ODEs, as part of a preconditioner, I get the system

$(A + D(t))x = b(t)$

where $A$ is a sparse matrix and $D(t)$ is diagonal. I'm currently solving this by taking the LU-decomposition of the system for each $t$. $A$ is constant throughout the solution.

Is there an algorithm that would allow me to solve this efficiently without recomputing a decomposition at every step?

  • $\begingroup$ I see some possibilities, but more information is needed from you. How accurately do you need to compute $x$? Are you using GMRES? Does $D$ and $b$ change slowly over time? What properties does $A$ have? What can you tell us about the original problem. Depending on your answers it may be fruitful to recycle preconditioners and/or Krylov subspaces and extrapolate initial guesses. $\endgroup$ Mar 23 '19 at 23:06
  • $\begingroup$ @CarlChristian OP writes that they are solving this system "as part of a preconditioner", so I assume they are interested in direct methods. $\endgroup$ Mar 24 '19 at 8:31
  • 1
    $\begingroup$ Ideally, I would be interested in a direct method; however, I can tolerate some inaccuracy in the solution. $A$ is the first order terms of a mass-action network. It has some very large off-diagonal terms. $D(t)$ captures second order terms of the mass-action network which have a special form for the problem I'm working on. $D(t)$ and $b(t)$ can change quickly. I'm using GMRES as part of the SUNDIALS suite. $\endgroup$
    – Sevenless
    Mar 25 '19 at 14:21
  • 1
    $\begingroup$ This is actually a frequently asked question: scicomp.stackexchange.com/questions/503 scicomp.stackexchange.com/questions/5101 scicomp.stackexchange.com/questions/2917 scicomp.stackexchange.com/questions/10630 scicomp.stackexchange.com/questions/2400 and the answer appears to be that there is no known way to use the fixed+diagonal structure. $\endgroup$
    – Kirill
    Mar 25 '19 at 21:20
  • $\begingroup$ @Kirill. I believe that you are trying help by strengthening the case against the fixed + diagonal structure. However, as this is OP's very first question here, I think we should mention that of the questions listed are all at least 5 years old, are not in the context of an ODE and do not include a time-dependence. It is precisely these differences which provide some options. $\endgroup$ Mar 25 '19 at 22:50

I do not believe that a direct solve is possible in your situation, there are simply too many applications with too many users which would benefit too much if it was possible. In particular, the problem of computing, say, the transfer function $T(s) = C(A - sI)^{-1}B+D$ of a linear time invariant system $$ \begin{align} x'(t) &= Ax(t) + Bu(t) \\ y (t) &= Cx(t) + Du(t) \end{align}$$ or solving the Lyapunov matrix equation $$AX + XA^T = BB^T$$ using the alternating direction implicit (ADI) method would be greatly facilitated by any algorithm which would allows to pass cheaply from a factorization of $A$ to a factorization of $A - sI$.

On the off chance that $D(t)$ has low rank, then the Sherman-Morrison-Woodbury formula would allow you to quickly solve systems with coefficient matrix $A + D(t)$ provided you can quickly solve system with coefficient matrix $A$. However, I gather that this is not likely.

I suggest that you view the solution $x$ of your linear system $$(A + D(t)) x = b(t)$$ as a function of $t$. It is likely that $x$ can be differentiate a few times with respect to $t$ which suggests that future values of $x$ can be extrapolated from past values. In particular, if $h$ is small enough, then $$ x_0 = x(t), \quad \text{or}\quad x_0 = 2x(t) - x(t-h)$$ represent simple and easy to implement initial guesses for the value of $x(t+h)$. More sophisticated approximation are possible, but they are not necessarily worth the effort.

I suggest further that your solve your linear system using a preconditioned Krylov subspace method. Instead of computing a new preconditioner for every timestep, I suggest that you recycle the preconditioner a few times and compute a new preconditioner only when the current one appears ineffective, i.e., the iteration count becomes too large.

I cannot promise that these options will give you the speed you are looking for, but at least they provide some parameters which can be tuned: the choice of the initial guess, the preconditioner and the number of time-steps before you compute a new preconditioner.


No (unless $D$ has few nonzeros). Variants of this question get asked a lot here, and the answer is invariably the same.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.