Solving system

This page deals with basics of solving square linear systems. For solving overdetermined and underdetermined systems, see Solving overdetermined systems and Solving underdetermined systems.

Basics

Solving a system $Ax = b$, where $A$ is a $n\times n$ matrix of scalars (complex, real, or otherwise) is possible if $A$ is invertible, and the solution is given theoretically by $x = A^{-1}b$.

Let $\hat{x}$ be the solution of $\hat{A} \hat{x} = \hat{b}$, where the matrix $A$ is disturbed as $\hat{A} = A + \Delta{A}$ and the right hand side as $\hat{b} = b + \Delta b$.

The absolute error of $x$ is $\Delta x = \|\hat{x} - x\|$ and the relative errors are defined as $\delta A = \|\Delta A\|/\|A\|$ and $\delta b = \|\Delta b\|/\|b\|$ and $\delta x = \|\Delta x\|/\|x\|$. We assume that relative errors in inputs are reasonably small.

If only $b$ is disturbed, the error is given by $$\delta x \leq \kappa(A) \delta b,$$ where $\kappa(A) = \|A\|\|A^{-1}\|$ is the condition number of matrix $A$.

In general, the error is given by

$$\delta x \leq \frac{\kappa(A)}{1 - \kappa(A) \delta A} (\delta A + \delta b).$$

Therefore, if the condition number $\kappa(A)$ is large, the accuracy of the solution of the system can deteriorate.

If $A$ is not invertible, Moore-Penrose pseudoinverse can be used to compute a "solution": $x = A^+b$.

Decompositions

Solving linear systems is not done by computing inverses of $A$ but by suitably decomposing $A$. Many decompositions are available, with $LU$ decomposition with partial pivoting being the go-to option. with $\frac23 n^3$ cost. Additionally, we have $QR$ and $SVD$ decompositions if more robust decompositions are desired, at a greater cost.

See full the list of Eigen's available decompositions: https://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html

If many systems $Ax = b_i$ need to be saved with only different right hand sides, the decomposition of $A$ can be computed in $O(n^3)$ and stored with subsequent solutions only taking $O(n^2)$ time.

If the matrix is symmetric and positive definite, Cholesky (called $VV^\T$ or $LL^\T$) decomposition can be used, with $\frac13 n^3$ cost.