# Speeding up solving \`Ax=b\` repeated when \`A\` is partially updated

**URL:** <https://discourse.julialang.org/t/speeding-up-solving-ax-b-repeated-when-a-is-partially-updated/92194>\
**Category:** Performance\
**Tags:** linearalgebra\
**Created:** [December 27, 2022, 6:00pm UTC](https://discourse.julialang.org/t/speeding-up-solving-ax-b-repeated-when-a-is-partially-updated/92194 "2022-12-27T18:00:45Z")\
**Posts on this page:** 1\
**Showing post:** 2

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [December 29, 2022, 11:52am UTC](https://discourse.julialang.org/t/speeding-up-solving-ax-b-repeated-when-a-is-partially-updated/92194/2 "2022-12-29T11:52:38Z")

</div>

> [@yhchang96](#):
>
> Are there better ways to solve for `x`?

There may be a better way, if the changes in `Amat` from one iteration to the next are “small” (in a sense made more precise below). Back in the stone age of computing (the 1990s), solving linear systems of equations was still pretty expensive even for matrices of order as small as a few hundred to several thousand. My fellow cavemen and I needed to solve many similar systems as we swept a parameter, frequency, over its analysis bandwidth. We adapted a well-known scheme called [iterative refinement](https://en.wikipedia.org/wiki/Iterative_refinement) for this purpose. Here is the basic idea of our modified iterative refinement (MIR) algorithm.

Suppose one wishes to solve the matrix equation Ax = b and that we have already computed the inverse (in practice, we’ll solve using the LU decomposition instead of explicitly forming the inverse) of a matrix A\_0 which is “nearly equal” to A. We compute a series of approximations to x as follows:

\begin{array}{ll} x^{(1)} = A\_0^{-1} b, & r^{(1)} = b - A x^{(1)} \\ x^{(2)} = x^{(1)} + A\_0^{-1} r^{(1)}, & r^{(2)} = b - A x^{(2)} \\ \hphantom{x^{(1)}} \vdots & \hphantom{x^{(1)}} \vdots \\ x^{(i)} = x^{(i-1)} + A\_0^{-1} r^{(i-1)}, & r^{(i)} = b - A x^{(i)} \\ \hphantom{x^{(1)}} \vdots & \hphantom{x^{(1)}} \vdots \end{array}

It is not hard to show that the i th residual r^{(i)} = \left(I - AA\_0^{-1} \right)^{i-1} r^{(1)}  
so that r^{(i)} \rightarrow 0 if the [spectral radius](https://en.wikipedia.org/wiki/Spectral_radius) of  
I - AA\_0^{-1} is less than 1. Inverting a matrix of order m (or obtaining the LU decomposition) is an O(m^3) operation, while solving the factored system and matrix-vector multiplication are only O(m^2). Therefore, this algorithm can save execution time if convergence is rapid enough (fewer than approximately m iterations are required).

Here is pseudocode for the algorithm from our [paper](https://ieeexplore.ieee.org/document/690832) on the subject.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/f/1fc9e7e3369769be4c7e38613e879fb4e9b24e9a.png)

For our application, this technique worked very well, resulting in substantial savings in execution time. Of course, YMMV. Hope it helps.

---

_[View the full topic](https://discourse.julialang.org/t/speeding-up-solving-ax-b-repeated-when-a-is-partially-updated/92194)._