Nonlinear system of equations whereas most of the equations are linear. How to minimise operation?

Question

Let us say we have a n * n system of equations like KU=F where K is a n*n matrix and U and F are n*1 vectors. K and F are defined and the final goal is to find U values. K is a sparse banded matrix and some of its components depends on U components. This dependence makes the whole problem nonlinear. Since the system of equations is nonlinear, we should use trial-error approach. We are using a decomposition method such as “LU” or “Cholesky” method to find U at each step. My question is about a case where 95% of the equations are linear which means that 95% lines and columns of the K are absolutely constants and only 5% of the equations are depending on U variables. Is there any way to avoid duplication of finding the LU or Cholesky matrices when most of matrix K is constant?

For instance, if we can find a part for the Cholesky matrix which is constant then we can start from that, and then only do the rest of the necessary operations for the 5% changing part of the equations. This will lead in to a lot of operation and time saving.

Any comment, and feedback, reference, ect is highly appreciated!

Best regards Meisam Jalalvand, PhD University of Bristol

Answer 1

Let us write your matrix equation as the following block format, $$\begin{bmatrix}A & B\\ C & D(x,y) \end{bmatrix}\begin{bmatrix}x\\ y \end{bmatrix}=\begin{bmatrix}b\\ c \end{bmatrix},$$ in which we have reordered the nonlinear portion of the problem into the block $D$. We highlight its nonlinearity by giving it arguments. As you previously explained, 95% of the rows and columns in your system are linear, so $D$ has very few rows and columns.

Considering doing block LU factorization, which results

$$\begin{bmatrix}I & 0\\ CA^{-1} & I \end{bmatrix}\begin{bmatrix}A & 0\\ 0 & S(x,y) \end{bmatrix}\begin{bmatrix}I & A^{-1}B\\ 0 & I \end{bmatrix}\begin{bmatrix}x\\ y \end{bmatrix}=\begin{bmatrix}b\\ c \end{bmatrix},$$ in which the nonlinear Schur complement $S$ is defined

$$S(x,y)=D(x,y)-CA^{-1}B.$$

We have block triangular matrices and block diagonal matrices, which can all be easily inverted to yield a direct solution,

$$\begin{bmatrix}x\\ y \end{bmatrix}=\begin{bmatrix}I & -A^{-1}B\\ 0 & I \end{bmatrix}\begin{bmatrix}A^{-1} & 0\\ 0 & S^{-1}(x,y) \end{bmatrix}\begin{bmatrix}I & 0\\ -CA^{-1} & I \end{bmatrix}\begin{bmatrix}b\\ c \end{bmatrix}.$$ Since $A$ is linear, its Cholesky or LU factorization can be precomputed, making matrix-vector products with $A^{-1}$ very easy. So every part of the decomposition can be precomputed except the inverse of the Schur complement, $S$. The linear part of $S$, namely the portion $CA^{-1}B$ does not change, and can also be precomputed. That final perturbation by $D(x,y)$ cannot be precomputed, but as you note, the block is small, so the cost of doing so is significantly less than redoing the entire LU factorization.

Going through the procedure requires ~5% effort for each new solution of the nonlinear equations, compared to the initial computation of the LU factorization for the subblock $A$, as desired.

Let me know if you have questions and I can refine this answer.