Using Matrix Differential Equations for Solving Systems of Linear Algebraic Equations

Abstract

Various ordinary differential equations of the first order have recently been used by the author for the solution of general, large linear systems of algebraic equations. Exact solutions were derived in terms of a new kind of infinite series of matrices which are truncated and applied repeatedly to approximate the solution. In these matrix series, each new term is obtained from the preceding one by multiplication with a matrix which becomes better and better conditioned tending to the identity matrix. Obviously, this helps the numerical computations. For a more efficient computation of approximate solutions of the algebraic systems, we consider new differential equations which are solved by simple techniques of numerical integration. The solution procedure allows to easily control and monitor the magnitude of the residual vector at each step of integration. A related iterative method is also proposed. The solution methods are flexible, permitting various intervening parameters to be changed whenever necessary in order to increase their efficiency. Efficient computation of a rough approximation of the solution, applicable even to poorly conditioned systems, is also performed based on the alternate application of two different types of minimization of associated functionals. A smaller amount of computation is needed to obtain an approximate solution of large linear systems as compared to existing methods.

1. Introduction

Exact analytic expressions in the form of infinite series of matrices for the solution of linear systems of algebraic equations were derived in [1] by integrating associated ordinary differential equations. These differential equations were obtained using a quadratic functional related to the system of algebraic equations and describe the orthogonal trajectories of the family of hypersurfaces representing the functional. More convergent matrix series were presented in [2] which can be applied to approximate the solution of the system of equations. Solution of linear systems based on the numerical integration of differential equations has originally been formulated in [3].

In Section 2 of the present book chapter, we use recently derived highly convergent series formulae for matrix exponentials [3] in order to construct improved iterative methods for solving approximately large systems of algebraic equations. In Section 3, we use novel functionals that allow to formulate differential equations which lead to a substantial increase in the efficiency of the solution process [4]. Independently of the starting point, the paths of integration of these equations converge all to the solution point of the system considered. At each step of the numerical solution, one passes, in fact, from one path to a slightly different one due to computation errors. The procedure does not require to find accurately an entire path but only the solution point which is common to all the paths. This is why we apply the simple Euler method [5] to integrate the differential equations. The computation errors are now determined by the magnitude of the second derivative of the position vector with respect to the parameter defining the location along the path. A related iterative method [6] is also described. In Section 4, two different kinds of minimization of the system functionals are applied alternately for quick computation of a rough solution of large linear systems [7].

2. Matrix series formulae for the solution of linear algebraic systems

Consider a system of equations written in matrix form as

$$\mathit{Ax}-b=0$$E1

where $A\in R^{n\times n}$ is an arbitrary nonsingular matrix, $b\in R^n$ is a given n-dimensional vector and $x={\left(x_1,x_2,\dots ,x_n\right)}^T$ is the unknown n-dimensional vector, with T indicating the transpose. Assume that $x$ is a continuous function of the real variable v over a certain interval and associate to (1) the vector differential equation of the first order

$$\frac{\mathit{dx}}{\mathit{dv}}=f\left(v\right)\left(\mathit{Ax}-b\right)$$E2

where $f\left(v\right)$ is a continuous function to be chosen.

2.1 Exact analytic expressions for the solution of (2)

Imposing the condition $x\left(v_o\right)=x_o,$ $x_o$ being a chosen position vector, (2) has a unique solution over a specified interval [4], namely,

$$x\left(v\right)=x_o+e^{-\mathit{Ag}\left(v_o\right)}\sum _{k=0}^{\infty }\frac{A^k}{\left(k+1\right)!}\left[{\left(g\left(v\right)\right)}^{k+1}-{\left(g\left(v_o\right)\right)}^{k+1}\right]\left({\mathit{Ax}}_o-b\right)$$E3

where $g\left(v\right)\equiv \int f\left(v\right)\mathit{dv}$ is a primitive of $f\left(v\right),$ i.e., $f\left(v\right)=\mathit{dg}/\mathit{dv}.$If $f\left(v\right)$ is taken to be $f\left(v\right)=1/v,$ then $g\left(v\right)=lnv.$ Choosing $v_o=1$ gives $g\left(v_o\right)=0$ and (3) becomes now

$$x\left(v\right)=x_o+\left(lnv\right)\sum _{k=0}^{\infty }\frac{{\left(Alnv\right)}^k}{\left(k+1\right)!}\left({\mathit{Ax}}_o-b\right)$$E4

Over the interval $v\in \left(0,2\right),x\left(v\right)$ can also be expressed in the form [3]

$$x\left(v\right)=x_o-\left(1-v\right)\left[I+\sum _{k=1}^{\infty }\frac{{\left(1-v\right)}^k}{k+1}\left(I-A\right)\left(I-\frac{A}{2}\right)\cdots \left(I-\frac{A}{k}\right)\right]\left(A{x}_o-b\right)$$E5

I denoting the identity matrix of order n. The solution of (1) is theoretically obtained for $v=0$, the series being in general extremely low convergent. Since the rate of convergence of the series in (5) is very small for v very close to the value corresponding to the solution of (1), i.e., $v=0$, this formula should be applied repeatedly for a v not too close to zero, $v\in$ (0,1), until a satisfactory small distance to the solution $x\left(0\right)$ is reached.

2.2 Highly convergent series iteration formulae

Practical formulae for an approximate solution of (1) can be derived by using matrix series that are much more convergent than the one in (5). This can be done by writing (4) in the form

$$x\left(v\right)=x_o+A^{-1}\left(e^{Alnv}-I\right)\left(A{x}_o-b\right)$$E6

and by expressing the matrix exponential in terms of series of superior convergence given recently in ref. [3]. Very close to the solution $v=0$, say for $v=e^{-N},N\gg 1$ , we have $e^{Alnv}=e^{-\mathit{NA}}$ and with [3]

$$\begin{array}{c}{e}^{-\mathit{NA}}-I=-{10}^{-q}{c}_q\mathit{NA}\left[I+\sum _{k=1}^{\infty }\frac{{\left(-1\right)}^k{10}^{-\mathit{qk}}}{k+1}\left(I+c_q\mathit{NA}\right)\right.\\ \begin{array}{c}\\ \\ \end{array}\times \left(I+\frac{c_q\mathit{NA}}{2}\right)\cdots \left(I+\frac{c_q\mathit{NA}}{k}\right)]\end{array}$$E7

where q > 0 and $c_q\equiv 1/ln\left(1+{10}^{-q}\right)$ , we get

$$\begin{array}{c}x\left(e^{-N}\right)=x_o-{10}^{-q}{c}_{q}N\left[I+\sum _{k=1}^{\infty }\frac{{\left(-1\right)}^k{10}^{-\mathit{qk}}}{k+1}\left(I+c_q\mathit{NA}\right)\right.\\ \begin{array}{c}\\ \\ \end{array}\times \left(I+\frac{c_q\mathit{NA}}{2}\right)\cdots \left(I+\frac{c_q\mathit{NA}}{k}\right)]\left({\mathit{Ax}}_o-b\right)\end{array}$$E8

In order to perform numerical computations, the number of terms in the series has to be appropriately chosen. To have a small $c_q$ in (8) one has to take a small q. For q = 0.5, e.g., $c_q=3.639409$. One may start with N = 15 which would require to retain about 50 terms in this alternating series to get a rough approximation of $x\left(e^{-15}\right)$. The computation is repeated with the new $x_o$ taken to be the preceding $x\left(e^{-N}\right)$ until an acceptable accuracy of the solution of (1) is reached.

Much more convergent formulae are constructed if we apply in (6) the expansion [3]

$$\begin{array}{c}{e}^{-\mathit{NA}}={\left(1+{10}^{-q}\right)}^{-p}\left\{I+p!\right.\sum _{k=1}^p\frac{{10}^{-\mathit{qk}}}{k!\left(p-k\right)!}\left(I-\frac{c_q\mathit{NA}}{p}\right)\left(I-\frac{c_q\mathit{NA}}{p-1}\right)\cdots \left(I-\frac{c_q\mathit{NA}}{p-k+1}\right)\\ -{10}^{-\left(p+1\right)q}{c}_q\mathit{NA}\left(I-c_q\mathit{NA}\right)\left(I-\frac{c_q\mathit{NA}}{2}\right)\cdots \left(I-\frac{c_q\mathit{NA}}{p}\right)\left[\frac{1}{p+1}I\right.\\ \left.\left.+p!\sum _{k=1}^{\infty }\frac{{\left(-1\right)}^k{10}^{-\mathit{qk}}}{\left(k+1\right)\left(k+2\right)\cdots \left(k+p+1\right)}\left(I+c_q\mathit{NA}\right)\left(I+\frac{c_q\mathit{NA}}{2}\right)\cdots \left(I+\frac{c_q\mathit{NA}}{k}\right)\right]\right\},\end{array}$$E9

$$p=1,2,\dots$$

The multiplication with $A^{-1}$ from (6) is avoided by arranging the first summation in (9) in terms of powers of A and by taking into account that

$$p!\sum _{k=1}^p\frac{{10}^{-\mathit{qk}}}{k!\left(p-k\right)!}={\left(1+{10}^{-q}\right)}^p-1$$E10

It is obvious that using (9) leads to computations with a much smaller number of terms retained in the series for a given N. With the same amount of computation we obtain an $x\left(e^{-N}\right)$ which is much closer to the exact solution $x\left(0\right)$ of the system of equations (1).

3. Methods of numerical integration

In this Section, we use special kinds of functionals that lead to the construction of differential equations which allow a substantial increase in the efficiency of the solution of large linear algebraic systems.

3.1 Vector differential equations and their application to the solution of (1)

Consider a functional of the form

$$F\left(x\right)={‖\mathit{Ax}-b‖}^{\alpha }$$E11

associated to (1) where $\alpha$ is a positive real number to be chosen. A real variable $v,v>0,v\in \left(v_o,v_S\right)$, is now defined by

$$F\left(x\left(v\right)\right)=F\left(x\left(v_o\right)\right)h\left(v\right)$$E12

where $h\left(v\right)$ is an appropriately selected real function with definite first and second derivatives in $\left(v_o,v_S\right),v_o$ corresponding to a starting point $x\left(v_o\right)\equiv x^{\left(0\right)},$ with $h\left(v_o\right)=1$, and $v_S$ corresponding to the solution vector of (1) $x\left(v_S\right)\equiv x_S,$ with $h\left(v_S\right)=0$. Denoting $F\left(x^{\left(0\right)}\right)\equiv F_o$, we have

$$F_o\frac{\mathit{dh}}{\mathit{dv}}=\alpha {‖\mathit{Ax}-b‖}^{\alpha -2}{\left(A^T\left(\mathit{Ax}-b\right)\right)}^T\frac{\mathit{dx}}{\mathit{dv}}$$E13

Thus,

$$\frac{\mathit{dx}}{\mathit{dv}}=\frac{F_o}{\alpha }\frac{\mathit{dh}}{\mathit{dv}}\frac{A^T\left(\mathit{Ax}-b\right)}{{‖\mathit{Ax}-b‖}^{\alpha -2}{‖A^T\left(\mathit{Ax}-b\right)‖}^2}$$E14

which is the differential equation to be integrated from $v=v_o$ to $v=v_S$. The second derivative of $x$ is

$$\begin{array}{l}\frac{d^{2}x}{{\mathit{dv}}^2}=\frac{F_o}{\alpha }\frac{d^{2}h}{{\mathit{dv}}^2}\frac{A^T\left(\mathit{Ax}-b\right)}{{‖\mathit{Ax}-b‖}^{\alpha -2}{‖A^T\left(\mathit{Ax}-b\right)‖}^2}\\ +{\left(\frac{F_o}{\alpha }\right)}^2{\left(\frac{\mathit{dh}}{\mathit{dv}}\right)}^2\frac{1}{{‖\mathit{Ax}-b‖}^{2\left(\alpha -1\right)}{‖A^T\left(\mathit{Ax}-b\right)‖}^2}\\ \times \left\{\frac{{‖\mathit{Ax}-b‖}^2}{{‖A^T\left(\mathit{Ax}-b\right)‖}^2}\right.\left[A^{T}A-2\frac{{‖{\mathit{AA}}^T\left(\mathit{Ax}-b\right)‖}^2}{{‖A^T\left(\mathit{Ax}-b\right)‖}^2}I\right]+\left(2-\alpha \right)I\}{A}^T\left(\mathit{Ax}-b\right)\end{array}$$E15

From (11) and (12) we get a useful relationship,

$$‖\mathit{Ax}-b‖={\left(F_{o}h\left(v\right)\right)}^{1/\alpha }$$E16

that allows to simply monitor the magnitude of the residual vector of (1) during the computation process.

As explained in the Introduction we apply the Euler method for the solution of (14) and compute successively

$$x^{\left(i+1\right)}=x^{\left(i\right)}+\eta {\left(\frac{\mathit{dx}}{\mathit{dv}}\right)}_{x=x^{\left(i\right)}},i=0,1,2,\dots$$E17

where η is the step size. In the absence of any hint about a good starting point $x^{\left(0\right)}$ corresponding to $v=v_o$, we have used the point along the normal from the origin $x=0$ to the surface $F\left(x\right)=\mathit{\text{const}}$ in (11) which is the closest to the solution point $x_S$ [8], i.e.,

$$x^{\left(0\right)}=\frac{{‖b‖}^2}{{‖A^{T}b‖}^2}{A}^{T}b$$E18

The function $h\left(v\right)$ and the parameter α are chosen such that the first and the second derivatives of $x$ in (14) and (15) remain finite when $v\to v_S$, i.e., when the residual $‖\mathit{Ax}-b‖\to 0$, while the errors evaluated with the second derivative are kept reasonably small along the interval of integration $v\in \left(v_o,v_S\right)$. For each $h\left(v\right),\alpha$ is determined by imposing the condition that the second derivative in (15) tends to zero as $‖\mathit{Ax}-b‖\to 0$. To decide on the value of $\mathit{\alpha }$ for a given $h\left(v\right)$, we require to have a good rate of decrease of $‖\mathit{Ax}-b‖$ at the beginning of the computation process, for instance to have (see (11) and (12))

$$\frac{‖A{x}^{\left(1\right)}-b‖}{‖A{x}^{\left(0\right)}-b‖}={\left(h\left(\eta \right)\right)}^{1/\alpha }\cong 0.8$$E19

when $\eta =0.1$. Finally, to solve (1), we determine the actual step size to be used for the numerical solution such that the errors at the beginning of the computation process are small, say

$$\frac{{\eta }^2}{2}\frac{1}{‖x^{\left(1\right)}‖}{‖\frac{d^{2}x}{d{v}^2}‖}_{x=x^{\left(0\right)}}<0.01$$E20

If the magnitude of the residual vector, which is computed at each step, does not decrease anymore the computation is continued with a new cycle of integration by applying (17) with $x^{\left(0\right)}$ replaced by a new starting point, i.e.,

$$\begin{array}{c}{x}^{\left(\mathit{new}\right)}=x^{\left(\mathit{\text{last}}\right)}-\frac{{‖{\mathit{Ax}}^{\left(\mathit{\text{last}}\right)}-b‖}^2}{{‖A^T\left({\mathit{Ax}}^{\left(\mathit{\text{last}}\right)}-b\right)‖}^2}\\ \times A^T\left({\mathit{Ax}}^{\left(\mathit{\text{last}}\right)}-b\right)\end{array}$$E21

where $x^{\left(\mathit{\text{last}}\right)}$ is the position vector from the preceding step. $x^{\left(\mathit{new}\right)}$ is the closest point to $x_S$ along the normal to the surface $F\left(x\right)=\mathit{\text{const}}$taken at $x^{\left(\mathit{\text{last}}\right)}$[8]. Subsequent cycles of integration are performed in the same way until a satisfactory approximate solution of (1) is obtained. If the difference between $x^{\left(\mathit{new}\right)}$ and $x^{\left(\mathit{\text{last}}\right)}$ is insignificant we find the point along the normal to $F\left(x\right)=\mathit{\text{const}}$, taken at $x^{\left(\mathit{new}\right)}$, where F has a minimum and, then, apply (21) again. It should be noted that as one approaches the solution point the direction of the gradient $A^T\left(\mathit{Ax}-b\right)$ tends to become more and more perpendicular to the direction of $x_S-x$ and, thus, the residual $‖\mathit{Ax}-b‖$ and the relative error $‖x-x_S‖/‖x_S‖$ will not decrease any more as expected. This is why the computation has to be continued by opening a new cycle of integration. For systems with higher condition numbers, this happens more quickly.

Numerical experiments obtained using $h\left(v\right)=1-v$ with $\alpha =0.45$ and also using $h\left(v\right)={\left(1-v\right)}^2$ with $\alpha =0.9$ shows that only two up to five integration cycles with a step size of 0.1 are needed in order to get an accuracy of about 1% for the solution of systems with condition numbers up to 100.

3.2 A related iterative method

The basic idea of this method is to find, starting from a point $x^{\left(0\right)}$, a point $x^{\left(1\right)}$ along the gradient of a functional (11) associated with the general system (1), such that the magnitude of the new residual vector is an established fraction of its initial value, $‖{\mathit{Ax}}^{\left(1\right)}-b‖=\tau ‖{\mathit{Ax}}^{\left(0\right)}-b‖,\tau <1$. Instead of performing an integration as in Section 3.1, one proceeds iteratively, i.e.,

$$‖{\mathit{Ax}}^{\left(i+1\right)}-b‖=\tau ‖{\mathit{Ax}}^{\left(i\right)}-b‖,i=0,1,2,\dots$$E22

for each iteration the starting point being the point found in the preceding iteration, with τ maintained as tight as possible at the same value. To do this, we impose that $F\left(x\right)$ in (11) varies from $x^{\left(i\right)}$ to $x^{\left(i+1\right)}$ in the same way for each iteration, namely,

$$F\left(x\right)=F\left(x^{\left(i\right)}\right)h\left(v\right)$$E23

where $h\left(v\right)$ is now a real function of a real variable, monotone decreasing in the interval $v\in \left[v_o,v_o+\eta \right],v_o\ge 0,\eta >0,$ with definite first and second derivatives, and with $h\left(v_o\right)=1$ such that

$$F\left(x^{\left(i+1\right)}\right)=F\left(x^{\left(i\right)}\right)h\left(v_o+\eta \right),i=0,1,2,\dots$$E24

Then, from (11),

$$\frac{F\left(x^{\left(i+1\right)}\right)}{F\left(x^{\left(i\right)}\right)}=\frac{{‖{\mathit{Ax}}^{\left(i+1\right)}-b‖}^{\alpha }}{{‖{\mathit{Ax}}^{\left(i\right)}-b‖}^{\alpha }}={\tau }^{\alpha }$$E25

and

$$\tau ={\left(h\left(v_o+\eta \right)\right)}^{1/\alpha }$$E26

$x^{\left(i+1\right)}$ is computed as

$$x^{\left(i+1\right)}=x^{\left(i\right)}+\eta {\left(\frac{\mathit{dx}}{\mathit{dv}}\right)}_{v=v_o}$$E27

i.e., with (14) where $F_o$ is replaced with $F\left(x^{\left(i\right)}\right)$,

$$\begin{array}{c}\begin{array}{c}{x}^{\left(i+1\right)}=x^{\left(i\right)}+\frac{\eta }{\alpha }{\left(\frac{\mathit{dh}}{\mathit{dv}}\right)}_{v=v_o}\frac{{‖{\mathit{Ax}}^{\left(i\right)}-b‖}^2}{{‖A^T\left({\mathit{Ax}}^{\left(i\right)}-b\right)‖}^2}\end{array}\\ \times A^T\left({\mathit{Ax}}^{\left(i\right)}-b\right),i=0,1,2,\dots \end{array}$$E28

in which, taking into account (21), one has to have $\left|{\left(\eta /\alpha \right)\left(\mathit{dh}/\mathit{dv}\right)}_{v=v_o}\right|<1$. This expression corresponds to the first step in the numerical integration by Euler’s method of the differential Eq. (14), starting from $x^{\left(i\right)}$ with a step $\eta$, or to the first two terms of the Taylor series expansion of $x\left(v_o+\eta \right).$

The starting value $x^{\left(0\right)}$ and the function $h\left(v\right)$ are chosen in the same way as in Section 3.1. For selected ratios $\eta /\alpha$ the iteration cycle continues as long as the residual $‖{\mathit{Ax}}^{\left(i\right)}-b‖$ decreases at a proper rate. Theoretically, to make $‖{\mathit{Ax}}^{\left(i\right)}-b‖=\epsilon ‖{\mathit{Ax}}^{\left(0\right)}-b‖$ with $\epsilon \ll 1$ one would need to conduct $\left(ln\epsilon \right)/ln\tau$ iterations. Since computation errors are introduced at each iteration (as in the Euler method), the initially chosen value of τ cannot be maintained the same as the iterative process continues. An approximate solution of (1) is obtained at the end of the iteration cycle as before, applying (21). Subsequent iteration cycles are performed with the starting point in each cycle being the point determined in the preceding cycle.

Numerical results were generated using, as in the method presented in Section 3.1, $h\left(v\right)=1-v$ but with $\eta /\alpha =0.2$ and $h\left(v\right)={\left(1-v\right)}^2$ with $\eta /\alpha =0.1$. Now, in both cases, $v_o=0,\left(\eta /\alpha \right){\left(\mathit{dh}/\mathit{dv}\right)}_{v=0}=-0.2$ in (28), and $\tau \cong 0.8.$ A substantial increase in τ, approaching $\tau =1$, or oscillations of its value during the first few iterations show that the computation errors evaluated from (15) (with $F_o$ replaced by ${‖\mathit{Ax}-b‖}^{\alpha }$) are too big. In such a situation, one has to decrease the factor $\left(\eta /\alpha \right){\left(\mathit{dh}/\mathit{dv}\right)}_{v=0}$ in (28), i.e., to decrease $\eta /\alpha$ for a given $h\left(v\right)$, which leads to an increase of τ and a decrease of $\left(1/2\right){\eta }^2{‖d^{2}x/{\mathit{dv}}^2‖}_{v=v_o}$. For accuracy of 1% for the solution of (1), one needs a number of three to six short iteration cycles, with about eight iterations per cycle, for systems with condition numbers below 100. Of course, as in the previous method, an increased number of iteration cycles is required to reach the same accuracy for the solution of poorly conditioned systems.

Remarks. One can easily see that $x^{\left(i+1\right)}$ in (28) can be expressed in the form

$$x^{\left(i+1\right)}=x^{\left(0\right)}+\sum _{\ell =0}^i{x}_d^{\left(\ell \right)},i=0,1,2,\dots$$E29

where

$$\begin{array}{c}{x}_d^{\left(\ell \right)}\equiv x^{\left(\ell +1\right)}-x^{\left(\ell \right)}=\frac{\eta }{\alpha }{\left(\frac{\mathit{dh}}{\mathit{dv}}\right)}_{v=v_o}\frac{{‖{\mathit{Ax}}^{\left(\ell \right)}-b‖}^2}{{‖A^T\left({\mathit{Ax}}^{\left(\ell \right)}-b\right)‖}^2}\\ \\ \times A^T\left({\mathit{Ax}}^{\left(\ell \right)}-b\right)\end{array}$$E30

with the solution of (1) given by the infinite series

$$x_S=x^{\left(0\right)}+x_d^{\left(0\right)}+x_d^{\left(1\right)}+x_d^{\left(2\right)}+\cdots$$E31

$x_d^{\left(\ell \right)}$ in (30) can also be written as

$$x_d^{\left(\ell \right)}=\frac{\eta }{\alpha }{\left(\frac{\mathit{dh}}{\mathit{dv}}\right)}_{v=v_o}\frac{{‖b^{\left(\ell \right)}‖}^2}{{‖A^T{b}^{\left(\ell \right)}‖}^2}{A}^T{b}^{\left(\ell \right)}$$E32

where

$$\begin{array}{l}{b}^{\left(\ell \right)}={\mathit{Ax}}_d^{\left(\ell -1\right)}-b^{\left(\ell -1\right)},\ell =1,2,\dots ;\\ b^{\left(0\right)}={\mathit{Ax}}^{\left(0\right)}-b\end{array}$$E33

This shows that the difference $x_d^{\left(\ell \right)}$ is a “rough approximation” to the solution of a system (1) whose right-hand side is, at each iteration, just the residual vector of the system in the preceding iteration, which decreases in magnitude from one iteration to the next. Thus, the method presented in this Section represents a practical, concrete implementation of the well-known idea of successive approximations/perturbations [9].

To search for a possible increase in efficiency, more general functionals of the form

$$F\left(x\right)=\mathcal{F}\left(‖\mathit{Ax}-b‖\right)$$E34

could be tested, where $\mathcal{F}$ and its first derivative are finite and continuous at all the points within the interval of integration. Now the corresponding (14) is

$$\begin{array}{c}\frac{\mathit{dx}}{\mathit{dv}}=F_o\frac{\mathit{dh}}{\mathit{dv}}{\left(\frac{d\mathcal{F}}{d‖\mathit{Ax}-b‖}\right)}^{-1}\frac{‖\mathit{Ax}-b‖}{{‖A^T\left(\mathit{Ax}-b\right)‖}^2}\\ \times A^T\left(\mathit{Ax}-b\right)\end{array}$$E35

Note. The paths of integration in the methods presented can be looked at as being the field lines of a Poissonian electrostatic field in a homogeneous region bounded by a surface of constant potential $F\left(x\right)={‖{\mathit{Ax}}_o-b‖}^{\alpha }$ and with a zero potential at the solution point, $F\left(x_S\right)=0$. In the particular case of $\alpha =2$ the ratio of the volume density of charge within the region to the material permittivity is constant, namely, $-2{\sum }_{i=1}^n{\sum }_{k=1}^n{a}_{\mathit{ik}}^2$ where $a_{\mathit{ik}}$ are the entries of A. By altering this electrostatic field one could eventually make quicker the approach to the solution point along the integration path.

4. Method of alternate minimizations

This simple method is based on the property of a functional of the form (34) associated with a general system of equations (1) to allow not only to minimize the value of the functional but also the distance to the solution point of (1). Using only minimizations along ordinary gradients of the functional is not efficient unless the system is very well-conditioned.

For computing efficiently an approximate solution of general, large linear systems of algebraic equations, we propose in this Section the alternate application of minimizations of a functional and of the distance to the solution point, along the direction of the gradient of the functional.

Consider the functional in (34) where $\mathcal{F}$ is a real function defined for all $x$ from a chosen starting point $x^{\left(s\right)}$ to the solution point $x_S$ of (1), monotone decreasing with $‖\mathit{Ax}-b‖,\mathcal{F}\left(0\right)=F\left(x_S\right)$. The gradient of $F\left(x\right)$ is

$$\nabla F\left(x\right)=\frac{d\mathcal{F}\left(‖\mathit{Ax}-b‖\right)}{d‖\mathit{Ax}-b‖}\frac{A^T\left(\mathit{Ax}-b\right)}{‖\mathit{Ax}-b‖}$$E36

i.e., in the direction of the vector $A^T\left(\mathit{Ax}-b\right)$.

The minimum of $F\left(x\right)$ along a straight line through $x^{\left(s\right)}$ in the direction defined by an n-dimensional vector d is found from the condition that $x-x^{\left(s\right)}=\mathit{\lambda d}$, where $\lambda$is a scalar to be determined, and $A^T\left(\mathit{Ax}-b\right)$ are perpendicular. This gives $\lambda$and $x$ for the minimum of $F\left(x\right)$, namely,

$$x_{minF}^{\left(d\right)}=x^{\left(s\right)}-\frac{\left[{\left(\mathit{Ad}\right)}^T\left({\mathit{Ax}}^{\left(s\right)}-b\right)\right]}{{‖\mathit{Ad}‖}^2}d$$E37

The point at which $F\left(x\right)$ is minimum along the normal taken at $x^{\left(s\right)}$ is determined by replacing d with $A^T\left({\mathit{Ax}}^{\left(s\right)}-b\right)$ in (37),

$$x_{minF}=x^{\left(s\right)}-\frac{{‖A^T\left({\mathit{Ax}}^{\left(s\right)}-b\right)‖}^2}{{‖{\mathit{AA}}^T\left({\mathit{Ax}}^{\left(s\right)}-b\right)‖}^2}{A}^T\left({\mathit{Ax}}^{\left(s\right)}-b\right)$$E38

The minimum distance between the solution point $x_S$ and a point $x$ along the direction of the gradient of $F\left(x\right)$ taken at $x^{\left(s\right)}$ is at

$$x=x^{\left(s\right)}+\mu A^T\left(\mathit{Ax}-b\right)$$E39

where the scalar μ is determined by requiring the distance $‖x-x_S‖$ to be minimum. This gives μ and $x$ for this minimum, namely,

$$x_{minD}=x^{\left(s\right)}-\frac{{‖{\mathit{Ax}}^{\left(s\right)}-b‖}^2}{{‖A^T\left({\mathit{Ax}}^{\left(s\right)}-b\right)‖}^2}{A}^T\left({\mathit{Ax}}^{\left(s\right)}-b\right)$$E40

which depends only on the residual vector and on $A^T\left(\mathit{Ax}-b\right)$ at the point $x^{\left(s\right)}$, independently of the form of $\mathcal{F}$.

As already mentioned, repeated minimizations using only (38) are not efficient for solving a general system (1). The same is true when using only (40). To obtain an approximate solution of (1), we apply the formulae (40) and (38) alternately, the starting point $x^{\left(s\right)}$ being each time the point determined in the preceding minimization. As in Section 3, when there is no indication about a convenient first starting point $x^{\left(s\right)}$, one can use the origin $x^{\left(s\right)}=0.$ Only a few iterations, up to ten, are needed for a solution with a relative error $‖\stackrel{\sim }{x}-x_S‖/‖x_S‖$ of about 1% for systems with condition numbers of up to 100.

The procedure is surprisingly efficient for a rough solution even for very ill-conditioned systems. For example, for the system $\mathit{Hx}=b$ where H is the Hilbert matrix of order eight and whose solution is $x_S={\left[1,1,\dots ,1\right]}^T$ [10], a solution of 6% accuracy is obtained with a starting point $x^{\left(s\right)}=0$ by performing only seven alternate minimizations (40), (38), with no equilibration or regularization preoperated on the system. By comparison [10], for the Gauss elimination method, the accuracy is only 40.6%, for the Gauss elimination with equilibration 9.15%, and for the Householder method with equilibration 5.6%.

Experimental results show that, as the new point $x$given by (40), (38) becomes closer to the solution point $x_S$, the direction $A^T\left(\mathit{Ax}-b\right)$ of the gradient in (36), i.e., of the correction terms in (40) and (38), becomes more and more orthogonal to the direction of $x-x_S$. This causes the magnitude $‖\mathit{Ax}-b‖$ of the residual vector of (1) and the relative difference $‖x-x_S‖/‖x_S‖$ to not significantly decrease anymore. To progress with the computation one can intercalate minimizations of $F\left(x\right)$ along directions that are different from that of $A^T\left(\mathit{Ax}-b\right)$. By numerical testing, it is observed that $x-x^{\left(0\right)}$, where $x^{\left(0\right)}$ is the original starting point, has at this stage, in most cases, a significant component along the direction of $x-x_S$. Thus, one can use such a minimization direction to try to improve the solution accuracy.

5. Conclusions

Highly convergent iteration formulae for solving general, large linear systems of algebraic equations are derived from exact analytic solutions of particular differential equations based on new, accurate series representations of the matrix exponential recently published in ref. [3]. Specialized differential equations which make it possible to monitor and control the computation errors and the decrease of the magnitude of the residual vector $‖\mathit{Ax}-b‖$ of (1) at each stage of the computation procedure are constructed and integrated numerically to approximate the solution of these systems. Two methods of the solution have been presented in this book chapter and the simplest Euler method is applied for the numerical integration of the vector differential equations. In the first method cycles of integration are used, each cycle starting from a convenient value of the unknown and continuing until the rate of decrease of $‖\mathit{Ax}-b‖$ becomes too small. The second method is an iterative method where a fixed rate of decrease of $‖\mathit{Ax}-b‖$ is imposed at the beginning of the iteration cycles. In this method, only the first step in the Euler method is computed at each iteration and each cycle of iteration is conducted until there is no significant change in the magnitude of the residual vector. These two methods are highly efficient for large systems with condition numbers below 100 since only up to six cycles with less than ten steps per cycle are necessary to get a solution accuracy of 1%, at each step within a cycle having to compute two matrix-vector multiplications. The method in Section 3.2 seems to be more efficient for systems with bigger condition numbers. The number of cycles of integration/iteration increases with the condition number and preconditioning should be done for ill-conditioned systems before attempting to apply the methods presented in this work.

The iterative method of alternate minimizations presented in Section 4 is intended for computing quickly a rough approximation of the solution of linear systems of equations. In this method, preequilibration or preregularization/preconditioning are not required to obtain useful results even for systems with poorly conditioned matrices.

The present Book Chapter has been intended to constitute a review of the work done so far on the subject matter. It describes the proposed new methods for an approximate solution of large linear algebraic systems using appropriately chosen matrix functionals and shows the procedures for constructing the concrete solution algorithms. At this stage, the results presented are only validated by preliminary numerical experiments which indicate the efficiency of the proposed procedures for deriving approximate/rough solutions of large systems. More numerical experiments involving systems with higher condition numbers, as well as theoretical results, will be presented in future work. It is my hope that other researchers will be attracted to this new area and rigorous theoretical results will also be established (theorems, etc.).