NATIONAL BUREAU OF STANDARDS REPORT

2076

SOLVING LINEAR ALGEBRAIC EQUATIONS

CAN BE INTERESTING

George E. Forsythe

<jS§>

U. S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS
U. S. DEPARTMENT OF COMMERCE

Charles Sawyer, Secretary

NATIONAL BUREAU OF STANDARDS

A. V. Astin, Director

cover of this report.

Electricity. Resistance Measurements. Inductance and Capacitance. Electrical Instruments. Magnetic Measurements. Applied Electricity. Electrochemistry.

Optics and Metrology. Photometry and Colorimetry. Optical Instruments. Photographic

Technology. Length. Gage.

Heat and Power. Temperature Measurements. Thermodynamics. Cryogenics. Engines and

Lubrication. Engine Fuels. Cryogenic Engineering.

Atomic and Radiation Physics. Spectroscopy. Radiometry. Mass Spectrometry. Solid State

Physics. Electron Physics. Atomic Physics. Neutron Measurements. Infrared Spectroscopy.

Nuclear Physics. Radioactivity. X-Rays. Betatron. Nucleonic Instrumentation. Radiological Equipment. Atomic Energy Commission Instruments Branch.

Chemistry. Organic Coatings. Surface Chemistry. Organic Chemistry. Analytical Chemistry.

Inorganic Chemistry. Electrodeposition. Gas Chemistry. Physical Chemistry. Thermochemistry. Spectrochemistry. Pure Substances.

Mechanics. Sound. Mechanical Instruments. Aerodynamics. Engineering Mechanics. Hydraulics. Mass. Capacity, Density, and Fluid Meters. Organic and Fibrous Materials. Rubber. Textiles. Paper. Leather. Testing and Specifications. Polymer Structure. Organic Plastics. Dental Research.

Metallurgy. Thermal Metallurgy. Chemical Metallurgy. Mechanical Metallurgy. Corrosion.

Mineral Products. Porcelain and Pottery. Glass. Refractories. Enameled Metals. Concreting Materials. Constitution and Microstructure. Chemistry of Mineral Products.

Building Technology. Structural Engineering. Fire Protection. Heating and Air Conditioning. Floor, Roof, and Wall Coverings. Codes and Specifications.

Applied Mathematics. Numerical Analysis. Computation. Statistical Engineering. Machine

Development.

Electronics. Engineering Electronics. Electron Tubes. Electronic Computers. Electronic

Instrumentation.

Radio Propagation. Upper Atmosphere Research. Ionospheric Research. Regular Propagation Services. Frequency Utilization Research. Tropospheric Propagation Research. High

Frequency Standards. Microwave Standards.

Ordnance Development. These three divisions are engaged in a broad program of research

Electromechanical Ordnance, and development in advanced ordnance. Activities include

Electronic Ordnance. basic and applied research, engineering, pilot production, field testing, and evaluation of a wide variety of ordnance materiel. Special skills and facilities of other

NBS divisions also contribute to this program. The activity is sponsored by the Department of

Defense.

Missile Development. Missile research and development: engineering, dynamics, intelligence,

instrumentation, evaluation. Combustion in jet engines. These activities are sponsored by the

Department of Defense.

® Office of Basic Instrumentation • Office of Weights and Measures.
NATIONAL BUREAU OF STANDARDS REPORT

ms PROJECT NBS REPORT

1101-10-!? 100 November 18, 19!? 2 2076

SOLVING LINEAR ALGEBRAIC EQUATIONS

CAN BE INTERESTING

by

George E, Forsythe

1

National Bureau of Standards

PREPRINT

1Sponsored in part by the Office of Naval Research, USN.

The publication, reprln

unless permission Is obi

25, D. C, Such permls

cally prepared if that <

Approved for public release by the

Director of the National Institute of

Standards and Technology (NIST) on October 9, 2015

part, Is prohibited

idards, Washington

irt has been speclfirt for Its own use.

CONTENTS

1* Introduction 1

2 . Survey of methods U

3. Linear iterative processes 9

h» Acceleration of stationary linear processes 1$

£ . Least-squares methods. 19

6. General descent methods 2U

7. Method of conjugate gradients 25>

8. Errors and "condition" 31

9 . Influence of computing equipment 314-

Bibliography 39
'

»

, -

„ ,

r » .

.

.

.

•
SOLVING LINEAR ALGEBRAIC EQUATIONS

CAN BE INTERESTING

1

George E. Forsythe

The subject of this talk is mathematically a

system of n linear algebraic equations in n

Ax - b

Here A is a square matrix of order n, whose elements are given real

numbers a^,. with a determinant d(A) / Oj x and b denote column vectors

and the components of b are given real numbers .(Complex numbers

would offer no essential difficulty.) It is desired to calculate the

"1 -1

components of the unique solution x e A b$ here A is the inverse

of A.

Such problems arise in the most diverse branches of science and

technology, either directly (e.g., the normal equations of the leastsquares adjustment of observations) or in an approximation to another

A talk delivered (under a similar title) before the Eugene meeting of the American Mathematical Society, June 21, 1952, by invitation

of the Committee to Select Hour Speakers for Western Sectional Meetings j received by the editors , 1952.

n Sponsored in part by the Office of Naval Research, USN.

1. Introdue^on.

lowly one. Consider a

unknowns, written

( 1)
.

.

.
2

problem (e.g., the difference-equation approximation to a self-adjoint

boundary-value problem for a partial differential equation) . These

two are very frequent sources of numerical systems 5 note that A > 0

(i.e., A is symmetric and positive definite) in both examples. The

order n is considered to range from perhaps 6 or 8 up to as large a

number as can be handled.^

^Steam [ 111], for instance, mentions the solution of a svstem of

order 2300 by the U.S. Coast and Geodetic Survey. The accuracy demanded

of an approximate solution \ varies widely; even the function which

is to measure the accuracy of \ varies or is unknown. Some "customers”

want to make the length |b - a|| small; some, | \ - A '4>|; others

have apparently thought only in terms of getting A"d exactly.

We all know that each component of the solution A ^b can be expressed as a quotient of determinants by Cramer's rule. We have all

evaluated determinants of orders 3 , U, and possibly £, with a^

integers; it is quite easy and rather boring. I therefore suspect

that the average mathematician damns the practical solution of (l)

as being both trivial and dull.

One defense of the numerical analyst (preface of [81]) is to

show that in many practical cases (say for decimal fractions a^ and

n » 10) getting A "H) efficiently (or at all) is actually not trivial,

but demands both know-how and planning. This point is especially

well taken against those whose program for solving (l) would be to

evaluate n + 1 determinants from their definition, employing nl multiplications per determinant. For n “ 26, for example,
(

.

'

'

I

, . ,

•

.

'

; ’

.

.

'

: . :

I

j
3

® Ofi

' (n + l) l * 10 , a number of multiplications which would take the

2 17 3

SWAC some 10 ' years. Actually, only about (l$)h: multiplications

2 , ^National Bureau 'of Standards Western Automatic Computer, an

electronic machine easily capable of 2000 multiplications per second.

See [6l] for a description, now out of date.

f

are needed to solve (1); see Bodewig [12], For n 26, (i/^n^ * 6000 .

and the multiplications would take the SWAC about 3 seconds.

It is my hope, on the other hand, to arouse the mathematician’s

t, interest by showing (section 2) the diversity of approaches to the

solution of (l), and by mentioning (sections 3 to 6) some problems

associated with selected iterative methods. The newest process on

the roster, the method of conjugate gradients, is outlined in section

7. Section 8 touches on the difficult general subject of errors and

"condition," while a few words are hazarded in section 9 about the

effect of machines on methods.

Whether or not the subject proves interesting, the bibliography

is intended to make the field look at least respectable l It is a

representative condensation of the 5>00-odd titles in the author’s

file, most of which are in [2y]. There are bibliographies on related subjects in Collatz [16] , Dwyer [20], Engineering Research

Associates [22], Frame [36], Frankel [37] s Franklin and Hankam [38a],

<vV AA<» Harvard Computation Laboratoiy [U7], Higgins [3&], KuroS, Markusevlc,

and Rasevskii [71], Motzkin [82], Schwerdtfeger [106], and

Taussky [116].
'
h

It is remarkable how little is really understood about most of

the methods for solving (l) . They are being used nevertheless, and

yield answers. This disparity between theory and practice appears to

be typical of the gulf between the science and the art of numerical

analysis.

2. Survey of methods. It usually surprises the uninitiated to

learn the variety of methods actually used to solve systems (l) . Surveys of methods are given by Bodewig [12], Dwyer [20], Faddeeva [23],

Forsythe [26] , Fox [33, 3k] , Frazer, Duncan and Collar [39], Hotelling

[ 58 ] , Householder [ 59, 60] , Zurmuhl [131], and others, but the spectrum

of known methods is considerably broader than any of these surveys

suggests. A classification of elimination methods is given by Jensen

[gy. A tentative classification of all known methods, with a bibliography of about U50 titles, is in [ 27 ] „

As mentioned in section 1, one can solve a system (l) explicitly

by use of determinants, while explicit solutions in other forms are

sometimes available^ Friedrich and Jenne [U0], for example, describe

the use of continued fractions.

The best known methods are based on systematic elimination of

unknowns from equations in the system (l) in the fashion of highschool algebra, as described by Gauss [l|2] . The elimination amounts

to triangularizing the matrix A by premultiplying it by a triangular

matrix, as Banachiewicz [ 5 ] and Turing [122] point out. The process

A-'”'—- ,—

can be rendered very efficient numerically by consolidating operations; see, for example, Benoit [7], Dwyer [20], and Turing [122].
'

*
When A is positive definite, the method is equivalent to the successive orthogonalization of the unit vectors in the A metric by the

Gram-Schmidt process [ 99]

With other metrics the orthogonalization

process yields different methods for solving (l) . All these elimination methods can also be performed on any matrix of submatrix blocks

formed by partitioning A, an idea suggested by Boltz [131 and based

on relations given by Schur [105>] . The various elimination methods

are direct , a term defined at the end of section 2.

There is a group of direct methods related to the characteristic

polynomial 0 of some matrix. For example, if one can learn that

0(A) s c n

An

+ **• + c^A +1=0, then, as Bingham [8] notes, one can

compute A~^b * -c

nA

n

"Hs - *°° - CgAb - c-jb. Similar remarks apply

when 0(H) is known for an operator H associated with the solution of

(l) j see section 1;. There are related methods involving the succes2 n

sive orthogonalization of the vectors Ax

q

, A x

q

, •••, Ax

q

in ihe

-1

I, A, A , or other metrics $ one of these is given in section 7.

There is an unpublished direct method of Hans Lewy using the

theory of congruences in n dimensions, applicable mainly when the

components of A~\> are integers. It is based on the use of stencils

for solving the system A x b , where a. . » 0 (if a. . is even) and

a^ “ 1 (if a^ is odd), and where b^ is similarly defined.

Of a quite different nature is a group of iterative processes

devised by Jacobi [63], Nekrasov [88], Richardson [98], and others,

and subsumed by Cesari [l£], Wittmeyer [12?], and Geiringer [UU] in

a general class of linear iterative processes to be discussed in
.

'
6

section 3. In these methods one starts with an arbitrary first vector

x

q

approximating A

-1

b. For k * 0, 1, 2, ••• , the components of the

k-th approximant vector x^, are systematically corrected, often in

the cyclic order 1, 2, ° °

- , n, but sometimes in blocks or otherwise.

After one cycle of the corrections has been completed, the components

of the resulting vector have beai obtained by solving a

matrix equation Bx^^ + Cx^ = b for x^-^, where B + C “A. The different linear processes are distinguished by the choice of B, which

must be an easily invertible matrix. For example, B can be the lower

triangle of A [88] , the diagonal of A [63], a sdalar matrix [98],

diagonal blocks of A (von Mises and Pollazek-Geiringer [123] , Hertwig

[U9]), etc. (Wien B is the lower triangle of A, the iterative method

is commonly called the ’Seidel process,' or the ’Gauss-Seidel process.’

But, as Ostrowski [91] points out, Seidel [10 7 ] mentions the process

but advocates not using it. Gauss nowhere mentions it. Nekrasov [88]

studies its properties and says it is Seidel’s process. It will

henceforth be called the ’cyclic single-step process.") Under appropriate conditions mentioned in section 3, x^. approaches A ^b, but

ordinarily no x^. equals A ^b.

Another group of iterative processes includes those of Gauss

[III], Seidel [IO7 ], Southwell [109, 110], Motzkin [2]^, and others.

These are called ’’relaxation methods" in [109, 110 ] . They are discussed by Black and Southwell [10], Fox [32], Temple [118] , and

others, and have proved especially useful in engineering work. They

are difficult to define precisely, since the computer is essentially advised to use all the art and artifice he can muster to
• - \ r ' * * ’ ; *
7

find x such that r « b - Ax is near 0. Their predominant feature,

however,, is that the components of x are corrected, not in a predetermined order, but in an order of "worst first If this feature

is adopted as defining relaxation, the iteration function (defined

in section 3) depends on x in a piecewise linear fashion, and the analytical character of the processes is completely different from that

of the related linear processes 0 Relaxation has been studied recently

in connection with solving systems of linear inequalities ; see

Agmon [2], and Motzkin and Schoenberg [83],

Other nonlinear iterative processes include the least-squares

iterative methods discussed in sections 5>,6, and 7. They start with

Cauchy [lU], and are synthesized by Temple [118], Rosser [unpublished],

Hestenes and Stein [53], and others* Special cases include certain

linear processes which essentially deal with a positive definite matrix.

For example. Kaczmarz [ 67 ] and Tompkins [120] interpret the system

(1) as restricting A ’Lb simultaneously to n hyperplanes . A first

guess x

q

is projected successively on each hyperplane in cyclic order.

Then the distance |x^ - A"° | decreases monotonically to 0. De la

Garza [19] proposes a related process. When A > 0, the cyclic

single-step method and all the relaxation methods are also leastsquares methods. This fact is the basis of many studies of the relaxation methodsj see [118] and Ostrowski [93].

Least-squares processes which are not linear include the gradient

,y methods of [llj.] developed in [118], and by Kantorovic [68] and Birman [9].
'

.

-

'
8

see section !?. Their culmination appears to be the conjugate gradient

method of Hestenes [£0, 5»U] , Lanczos [72] ? and Stiefel [llh, 5k], this

is a finite iteration described in section 7 . A gradient method' in

a more general metric is mentioned in section 6 .

Seemingly very different are the Monte 'Carlo methods for solving

(l), employing random sampling of a certain chance variable whose expectation is A \>. One such method, devised by von Neumann and

Ulam [unpublished], is described by Forsythe and Leibler [28 ]

5

Wasow

rig] devises another Monte Carlo method. Both are based on propert ties of discrete random walks in a space of n points. When the system

(l) represents Laplace's or related difference equations in one or

more dimensions, the Monte Carlo methods have a longer history; see

the exposition and bibliography in Curtiss [18]. The methods are

theoretically fascinating, but there is little evidence yet of their

practical utility for solving linear systans.

The iterative processes for solving (l) are likely to converge

slowly, and a number of tricks have been devised to speed them up,

called acceleration processes. Accelerations of linear processes may

themselves be either linear or nonlinear; some are described in

section U. A number of processes for accelerating the nonlinear gradient methods are mentioned at the end of section and the conjugate

gradient method of section 7 may be considered also as an acceleration.

1 The distinction between direct and iterative methods is ordinarily stressed in numerical analysis; see [23], for example. Applied

to systems (l), a direct method is one which yields k~\> exactly
' < « ; ' • . . n
9

after a finite number of arithmetical operations,, if the latter are

performed without round-off error. An iterative process can ordinarily

yield A~"ho only as the limit of a sequence of exact arithmetical operations. However,, it must be remembered that as soon as calculations

are rounded off (as ordinarj^r occurs in machine calculation), direct

methods disappear except for quite simple problems, and all methods

become iterative! see the end of section 8. The practical distinction

between methods for solving (l) then appears to depend on: (a) the

speed of the convergence to A~'H>, and (b) the simplicity of the computation at each stage of the iteration. A two-way classification of

methods might well be based on their behavior with respect to properties (a) and (b) , and that of [2£] was roughly of this type. Two

difficulties make such a classification imprecise. First, the theory

of round-off error is too poorly developed to yield rates of convergence to A”^b. Second, the practical criteria of simplicity in

machine computation vary too greatly among different machines.

One may also distinguish whether the matrix A is altered in the

course of solution, as in elimination, or whether it is retained in

its original form, as in the conjugate gradient process. This is

probably a crucial distinction in practice, for, when the original A

is retained, the round-off errors seem to accumulate more slowly!

see section 9.

3* Linear iterative processes. An iterative process for solving (l) (or other equation) is determined by the functions wherewith the (k + l)-th approximant to A*^b, x

k+l-’

is derived from the
<

'

-

msm
10

earlier approximants XQ?Xp •••

, x^.. If the only argument of F^ is

x, s the iterative process is said (Schroder [103]) to be of first

degree; “ F^(x^) . If the function F.^ is independent of k, the

process is called stationary . A bibliography on iteration as such has

been prepared by Schwerdtf eger [106]

„

As elsewhere in mathematics, the most studied functions are the

linear ones . We Introduce the (most general) linear iterative process

of the first degree by the definition

( 2) *k +!

" W Vic + \ >

where the H^. are square matrices and v^ are vectors. If the iterative

process described by (2) is to solve (l)

,

it seems essential that the

solution A ^b be a fixed point of F^:

( 3 ) A

J

'b “ F^(A ^b) « H^A + v^

We demand that (3) hold* It follows that

(U) \ (i - v a“^ * v •

If Hj and M, are independent of b, then (6) follows from (U) . Thus:

'RtJvx — The most general stationary linear iterative process of the

first degree for solving (l) which is independent of b and which

satisfies (3) Is defined by the functions

(*) Vi ' W Vic + V »

where the square matrices and depend only on A and satisfy the relation
'

’
11

(6) + \A - I *

If the process is stationary , then

(7) ^

+1

53 + Mb ,

where the square matrices H and M depend only on A and satisfy the relation

( 8 ) H + MA - 1 .

Stationary processes of the type ( 7 ), ( 8 ) have been studied by

Cesari [15], Wittmeyer [127], Geiringer [UU], and others. They include the cyclic single-step iteration of Nekrasov [ 88 ] and Liebmann

[7U], and those of Jacobi [63], Richardson [98], Frankel [38], Young

AA^ /W yw

[129], and many others. In the cyclic single-step process, for example, one writes A as the sum of two matrices B “ (b. .) and C (c. .),

where

r

a

ij (i 4

3)

(i <j)

and c. . ® ' 1J

o (i » j)

a

ij

^ < 3)

Then + CXy b, or + B~o; it is assumed that no

a.. - 0. Thus H - - B

_1

C, while M - B"1

.

Let us consider the convergence of the linear processes, under

the assumption that all arithmetic operations are carried out with

perfect accuracy. From (2) and (3) it follows that

(9) Vi - A_lb A_lb)

>
>

.

'

"
12

whence

(10) 2^

- iTh -‘Kk

(x

o

- Ch) ,

where K^. » ' ' * H

iH

0

* or<^er that - ATM^ 0 for

arbitrary x

q

it is therefore necessary and sufficient that

(11) lira K, z - 0 , for every vector z. k-»oo K

In practice condition (11) is usually known to hold only in

certain special cases, such as when?

(a) all are polynomials in some one matrix - for example, A;

(b) all = H (stationary process) |

(c) with some norm function ||*|| , all f| H.

|J 6 1 - £ <1.

Henceforward we consider only case (b)

stationary linear processes. Then » H^, and it is known (see [ 9U], for example) that (ll)

holds if and only if each eigenvalue X^(H) of H is less than one in

absolute value. Thus, we have derived the well known result that in

the stationary linear process x^ A“

1b for all x

q

if and only if

all | A

±

(H) | <1.

This result, while exceedingly important, hardly begins to solve

the practical computer’s problem. He raises the following questions:

(i) How can one tell whether all j X ^(H) j < 1 ? (ii) If all

| X^(H) J < 1, how fast will the convergence be? (iii) If the convergence is intolerably slow (as the computer realistically expects),

how may it be speeded up? To a considerable extent these questions
.
13

are answered by a knowledge of the asymptotic behavior of the error

vector - A ^b as k 4 oo„ Question (iii) will also be dealt with

in section U.

The general theory of the Jordan canonical form [ 78 , 9U] of the A/ '\Sm* matrix H can be applied to discuss the asymptotic behavior of

x

k

- A

1b. Suppose, as is usually the case, that H has the diagonal

canonical form

(12 ) H A-c.r? 111 s

T where c^, r^ are respectively the column and row eigenvectors of H

belonging to . Suppose that x

q - A H> * From (10) one

then finds that

(13) n

" A"^ - H

k

(x - A

-1

b) - J

i 5®!

If X n

is a unique eigenvalue of maximum modulus, and if y / 0, then

( 1U) x^ - A~^b - Yn X^c

n

(as k -4 00 )

This indicates that asymptotically i A

^

along a certain line, a

fact which is often useful in accelerating the convergence of

When An

is unique but « 0, the relation (lU) does not hold

for exact operations* However, a calculation involving round-off

errors will soon have the practical effect of making y^ / 0. If

several eigenvalues dominate in absolute value, formula (lU) must

usually be modified to include several terms . In any case, when (12)
'

*

*
iu

holds, formula (13) shows that |x^ ~ A^^b j approaches zero linearly —

i.e., with the order of the k-th power of rftax^| A^|. When (12) fails

(i.e., when H does not have a diagonal Jordan canonical form), formula

(13) has to be altered. The important fact, however, is that the

basic mathematical theory is already in existence for dealing with

linear iterative processes — at least as long as round-off errors

are not considered in detail. (It is commonly supposed that a nondiagonal canonical form would never arise in practical analysis. However, if one uses the Liebmann [?U] process to solve the usual difference equation corresponding to the Dirichlet problem in one or more

dimensions, the matrix H turns out to have a non-diagonal canonical

form.)

Recent developments in stationary linear processes have included

a suggestion of Aitken [U] and others for making H > 0 (i.e., a positive definite symmetric matrix), so that all A^(h) will be real and

positive, insuring that (13) will in fact be dominated by one term.

In the cyclic single-step process this can be achieved by solving for

the components in the " to-and-fro*' order 1, 2 9 ••• , n, n - 1, •••

, 2

and repeat.

1 According to Ostrowski [91], for the cyclic single-step prbcess it

was first proved by Pizzetti [95>] that, if A > 0, all i A ±

(H) ! < 1. A

converse recently proved by Reich [97], and more simply by Ostrowski

r\s~^ 1 [92], asserts that if A is symmetric and each a^.^ > 0, then all

v'VV* ““ —— “ H

| X^(H) I < 1 if and only if A > 0. Collatz [17] > Stein and Rosenberg [113], and others have studied the relations between the A

^ (H)
.

l' ; . * - ’ - ' % '
35

for the cyclic single-step process and those for the related total-step

process of Jacobi [63]. Other recent developments [38. 129] include

the alteration of H by simple devices so that max^j A is reduced as

much as possible. In the Liebmann process for the Dirichlet problem,

a certain systematic "over-relaxing" has just this effect.

U. Acceleration of ^tatioiw^- linear proc esses . Since ordinarily

the stationary linear process (?) seems to converge slowly, or even

diverge, it is of the greatest practical importance to devise methods

to improve the convergence. Any procedure which replaces by a

vector closer in some sense to A

=^b is loosely called an acceleration

of the iterative process ( 7) . An acceleration is ordinarily a

summability

best understood accelerations are, like the iterations themselves,

linear processes , in which the improved vector is a linear combination

of V x

l>

'**

> V

The simplest linear acceleration is useful when H has a unique

dominant eigenvalue A

n

? By (lU),for large k

applied to the sequence defined by (13) . The

05) Vi ' A_lb 1 X

n

(x

k

" A_lb) •

Hence the formula

(16)

-L . y

“k*i" * nSc

D a n

may be expected to give a better approximation to A "Hd than gives,

even when A

n

is known only approximately. This idea has been applied by Flanders and Shortley [2k], Lyusternik [77]* and others.
,

-

, . .

- .

.1 ,
16

The same idea can be extended to take care of two or more dominant t A.. Let E be the operator increasing by one the subscript k of x:

x^+^

* Ex^. One can write (16) in the form

(17) » P

1

(E) ^ ,

where P, (E) " (1 - A (E - A ) is a polynomial of degree one in X XI XX E • Let m be an integer . Consider the formula

(18) A

-1

b i P

m

(E) Xy ,

where P^ is some polynomial of degree m. How shall we determine P^

so that (18) is a good approximation? If we can answer this effectively, then (18) is a useful linear acceleration of the process (7).

In order that the operation (l8) not worsen an already good approxi

mation x^, we must have

(19) P(l) “ 1 . m

For if x

k

- pT\} then, by (3), x^+1 * x

k+2 -

••• - k~\, so that

VE) \ *k ' p«

(1) ch> -

If m - n the order of H, and if 0(A) is the characteristic

polynomial of H, then the choice P^(E) * 0(E) /0(l) is perfect, in

that P

m

(E) x^ is exactly equal to

0 is unknown, and to obtain it would involve precise knowledge of

all the! eigenvalues A^ of H.

Suppose that the eigenvalues although not known precisely,

are known all to be in a certain closed set R of the complex plane

for all x^. But ordinarily
. >.

'
17

(for example, a real interval). By (8), since d(A) ^ 0, the number 1

cannot be an eigenvalue of H. Hence xire may assume that ft is bounded

away from the point 1. Now the eigenvalues of P

m

(H) are

Since, by (13),

P

m

(E) ^

- A~\ - P

m

(H)(x^ - , ! 0

it is essential for the success of the acceleration that all |P (A.) I

1 m i1 be small. This leads to the approximation problem of determining the

polynomial P

m

of degree m such that (19) holds and such that

( 20) max |pj A ) I is a minimum

X£R “

Such polynomials, named generically after Cebysev (* Chebyshev *

Tschebyscheff * • • •) arise frequently in numerical analysis.

If R is a real interval not containing 1, the minimizing polynomial

P^(X) is proportional to the ordinary Cebysev polynomial T ( A ) for

the interval Rj see W, Markoff [79]. For just this reason the polynomials T^ have been used several times in matrix problems to suppress

the effects of unwanted eigenvalues; see Abramov [1], Flanders and

Shortley [ 2U3

,

Gavurin [b3], and Lanczos [73]* A recent review of

y' ^ Cebysev approximation problems in the complex plane is in [80].

The real difficulty in numerical analysis is that R is not knoxm.

How can the information gained from computing x

q

, x^, •••, x^., •••

, x^

+m

V y/ be used effectively to localize R, so that the Cebysev problem can

be formulated? After R is localized, how can the corresponding
.

.

-

'

' ,

;

„
18

P^(E) be determined approximately? These are important questions

which are in need of treatment. 'The use of symmetric matrices H helps

a great deal, by restricting R to the real axis.

Another acceleration device is due to Schulz [lOU], Suppose

A » I - H, and that (2l) converges, A linear process analogous to

(7) obtains A

^ as the sum of a Neumann series

(21) (I - hF1 - I + H + H

2 + •••

,

often slow to converge. By [10U] one can obtain the partial sum n

X =I + H+"'*+H2 1

of (21) by n iterations (22 ) s

( 22) X

k*l

’ X

k<

21 - V » X

o ' 1 •

a total of 2n matrix multiplications. The use of (22) to get A is

sometimes called Newton’s process. Can this idea be adapted to solving

the system (l) without getting A~^ first?

For accelerating the convergence of any linear process one also

has the S ^-process of Aitken [3$ U] and its extensions by Shanks [108]

rv/-. and Samuelson [101], Lubkin [ 76 ] has studied it as a nonlinear sequence-to-sequence summability process. It requires no knowledge of

the A . (H) . Let y. represent an arbitrary, but fixed,, component of

x^.. Then the functional character of y^ in a linear process is

given by

(23 ) n

lx 00 i=l

+ V 8 . A?

1 1

where y^ is the desired component of A "Sd, and s^ are numbers.
'

• . v

•

1

. .. .

'
19

To determine y from (23) it is theoretically sufficient to have exactly 2a + 1 successive values of y^ — for instance,

y , y , **•

, ^n 9 Prac^ce the elimination of the and A^

would be too tedious, but frequently one A^ ? say > , predominates

in (23) , and for moderately large k one may ignore s^, for

2 ' i “ 1, , n - 1. In the 8 -process one assumes that, with suffiV cient accuracy, y^. ® y^ + sn

A

n5

and calculates y ^

from three successive y^. The formula is thiss

¥w'Vi _

( V

2

y »

' 7k " 1^7

•

In [108], m exponential terms are eliminated from (23) at once.

If m « 2, for example, one can show that

f?. Least^squares methods . A variational approach to solving (l)

is very popular now. In a general treatment over the complex field

(the present treatment is confined to the real field), Hestenes and

Stein [£3] take a matrix R > 0 and let (2 Rz) 2 » |z| R

be the R-length

of a vector z. (The T stands for transposition.) If B " A RA and

c A^Rb, then B > 0 also. Let the deviation of x from A be

measured by
. . ? ' » ' - .. : fcfci-rtt-
20

(2U)

f(x) * | Ax - b|^ » |x - A

=1b|g

* x

TBx - 2xT

c + jb |j?

Starting from an initial vector x

0 , one can move x in various directions with the object <f minimizing f(x). Clearly the minimum is attained just when x - A~'

!'b.

To simplify the exposition we now assume A > 0 and take R * A

-

^.

Then

(25>) f(x) - x

TAx - 2xTb + b

T

A~

1b .

Although f(x) is not computable unless A "Ha is known, it is sufficient

in practice to minimize the computable function f(x) - b^A”''‘b. Fix x,

and let r » b - Ax be the residual of the systan (l) at x. Let d / 0

determine a direction. Since f(x + °< d) °<- d Ad - 2<xd r + f(x),

the value of the real parameter for which f(x + ^-d) is a minimum is

( 26) oCf - d^r/d~Ad

this is called the optimum oC (corresponding to x and d)

.

For any ©C » one can compute f(x + <vd) from the relation

f(x) - f(x + d) - (2«-

ot ^)d^Ad, i.e., f(x) - f(x + /3°<^d) -

(3(2— /S)«-’

r^dTAd. Thus (when d^r / 0), f(x + ^^^d) < f(x) for

0 < /3 < 2. The greatest reduction in f(x) comes when - 1.

i We now describe a general least-squares iterative process for

solving (l) . There must be prescribed? (i) a start x

o

$ (ii) a sequence of directions ; (iii) a sequence of ratios W-

(In but not here, R must also be prescribed.) For each
-

h

.
21

k a 0, 1^ 2, ***

j one determines °C ® by (26) so that f(x^. + ^ d^)

is minimized. Then one lets

(27) Vi m \

* 'V'fo >

T / T where r^ - b - Ax^. and “ - d^r^/d^Ad^. (It is not

the A and the

If the sequences

in [f>3] that, independently of x

q

, f(x^) -*-» f(A*”'

Lb),

that » A "H). The conditions are;

d^, be determined a priori; they may

and fdjT satisfy (28) and

specified that

depend on x^.)

(29) it is shown

as k —» co, so

(28) 0 < &

1 < /3

y

4 2 - ^ < 2 (all k) }

(2?) 0 < S

2 - |d^rk

|/(|dj

c |-!rk

|) (all k) .

In [E>3] an alternate to (28) states that the d, recurrently span the

space in a certain uniform manner.

Among others, the following two least-squares processes are also

linear iterative processes; (a) the and are independent of

the x^j (b) d^ r^ but &

* ^k°^k * a constant. Wien in

(a) all

/?k

® 1 and the d^. are the coordinate unit vectors in cyclic

order, one has the cyclic single-step process of Nekrasov [88], The

use of f(x) is veiy useful in studying any single-step process. In

(b) one can write the linear iteration function of section 3 in the

form F(x) * x + ol(b - Ax) = (I - °CA)x + b; the process is due to

Richardson [98] and to von Mises and Pollazek-Geiringer [123], and

converges whenever 0 < °^ < 2/max. ^ .(a).

ZL ZL
'

1
22

In the general case, however, the least-squares process is nonlinear. When 83 r^ 53 - \ grad, f(x^), it is called a gradient

method (or method of steepest descent) . When f3^ 3 1, one has the optimum

gradient method (since f(x^ + d^) is minimized as a function of ex),

proposed by Cauchy [lU] and studied by Temple [118], Kantorovic [68],

Birman [9] s and by Hestenes and Karush [5>2] for the eigenvalue problem.

Some variations of the method are treated by Krasnosel’skii and

Krein [70a].

By (2£) the surfaces f(x) constant are similar and similarly

situated ellipsoids whose common center A “4) we seek. Ary approximant x^. lies on a certain ellipsoid of the family. The gradient,

-2r^, lies in the normal to at x^, Now x^^ is the unique

point of (x^.) for which f(x) is minimized. Since f(x) is a quadratic function of distance along the normal, x^.+^

is located halfway

between the two intersections of '^C25^

“with S^. Moreover, x^

+ -^

i-s

the point where ITj (x^.) is tangent to an ellipsoid of the family.

Let 0 < A * * -

[68] shows that

*= be the eigenvalues of A. Kantorovic

(30) Fi < 1 (all x^)

From this it follows that f(x^) i 0, so that x^. A ^b and the process converges. As we mentioned in section 3 } the practical computer’s

problem merely begins with such knowledge. Some experience convinces

him that (a) is frequently very close to 1, and (b) after a few
_

•'

.

'

>

.

.
23

steps [f (x^.+^ )/f becomes and remains very close to see [25] .

Observation (a) is borne out by the remark [6^ p. 591 that, for cerm tain matrices A of type C

1 C, X

n/A ^

is likely to be of the order of

2 -2 n , whence 1 - ** n . For finite-difference approximations to

the Laplace operator over rectangular regions Frankel [38] shows that

^ , where n is the number of lattice points. As to (b) , 2 for any A the value can always be attained by ^(x

jc+^)/;8( xic

) when

\

^ assumes certain directions in the plane spanned by and

u

n

(defined below). From (a) and (b) we conclude that the optimum

gradient method is frequently too slow for practical use.

As a guide to the possible acceleration of the method it would

be valuable to know the asymptotic behavior of x^ - A~^b, if arithmetical operations are done exactly. But, because x^+^

is obtained

from by a rational cubic transformation, theorems are hard to

prove*. It is a conjecture of Forsythe and Motzkin [30], proved only

for n • 3 in [ 29 ], that in the optimum gradient method the error

vector x^. - A ^b is asymptotically a linear combination of the eigenvectors u

n, u^ of A belonging to the largest ( X and least ( X

eigenvalue of A. (If there are eigenvectors of A orthogonal to

x

q - A~^b, one disregards the corresponding eigenvalues in determining X and X

n

„) A proof of the conjecture for n U would be very

desirable because, when the conjectured asymptotic relationship

holds, for all sufficiently large k the points x^, x^-^, x

k+2 f A ^b are asymptotically coplanar. Thus one could accelerate the

convergence of the optimum gradient method by occasionally minimizing
.

1

*

'

"

.

.
2k

f(x) in the plane through the endpoints of x^., and x^

+ ^

e The

method has been used successfully' in experiments [25] for n * 6 on

•~\S\ an IBM Card-Programmed Calculator, where the average value of

f(Xk+i)/f (x^) over a series of about 100 steps was reduced from .9733

(optimum gradient method) to ,62[|5 (optimum gradient method with acceleration every ninth step) .

A second idea to speed up the optimum gradient method, the

ft -optimum method proposed by Hartree [U6] and by Hestenes and Stein

[£3], is to give ft, some constant value ft in the range 0 < ft < 1.

In the test problem of [25], Stein [112] finds ft “ 0.9 to be approxir\s\j

mately best, and with it the average value of f(x^. +^)/f (x^.) is found

[2^] to be .820U| in [112] .8065 is found for a shorter run. Although

/\ry/

the convergence of the ft -optimum method is slower (in this test) than

that of the accelerated optimum gradient method, the former has the

considerable advantage of being shorter to code for automatic machinery.

The success of the ft -optimum method is perhaps due to an inherent

instability of the transformations . In 2 dimensions the transforma-

-1

to ~ A b has no stable fixed directions when

ft is slightly less than 1 [Motzkin and the author, unpublished]

.

6. General descent methods. The function f of (2h) has the

tion of

*k

following significant properties: f(A ^b) 0| f(x) > 0 if x / A \>J

f is a convex function of x. Any f with these properties might serve

as the basis for what Professor Motzkin calls a descent method. In such

a method one has a first guess x

q

, and sequences {v} and {X} as

in section 5. As before, one finds a £ minimizing f(x^ + « d^)

,

and selects x^^ » x^. + /d^ ^d^

.
.

,

.

'
2?

Other suitable f would be

n

( 31)

(32) f(x) - max Jiv

i

Methods employing the latter norm functions f, with d^ “ - grad fCx^.),

are somewhat related to - but apparently do not include - the piecewise linear iterative processes. Agmon [2] discusses a relaxation

method for linear inequalities from somewhat this point of view.

Zuhovickii [130] gives a gradient method in the metric (32) for

minimizing f(x) for incompatible systems Ax b.

7. Method of conjugate gradients. It is easy to show that the

acceleration step of [25>], discussed in section 5, is equivalent to

finding the unique point x

f for which f(x') - f(x^ + Q

r^ + ^Ar^)

is minimized as a function of the two parameters <x qand This

V’ suggests a generalization to p parameters discussed by Kantorovic

[68] and Birman [9]. Let x be any point, and let r * b - Ax . Define f(x) by (25). Let if (x

Q

) be the p-space of points

x + c<r + d,Ar + ... + oc xr . where &•••&- _ are

o o o 1 o p—1 o

5 o’ 3 p—1

real parameters. Let x-, be the unique point in ft for which f(x)

is minimizeds

To determine x^ one has to solve a linear system of p equations in

I the p unknowns 9 ••• ,

., .
.

.

.

.

.

'
26

Define A ....

3 A as in section 5> • Let the interval

n A 4 A ^ A be transformed into the interval - 1 4 t 4 1 by the

j. n d affine transformation t ® Z ( A ), carrying into -1, and into 1.

Then Z (0) - S * ( X

n

+ A^)( A

r

- A^) Let Tn

(t) be the ordinary

Cebysev polynomial of degree n, normalized so that Tn

(l) max

]_4t6l

I

T

n

(t) I TM 1. In these notations Birman [9] has proved that, for each A,

l ~~

'2

' TTT7 ‘ Pp < 1 (al1 x

o>

•

i 1

*' (3b)

f(Xj)

Thus (3b) is the extension to p > 1 of (30)

.

Let t^ tT(X^), so that t^ “1, t

n

“ -1. For a given A, the

value

pp

cannot necessarily be attained by |?(x^)/f (x

Q )p . It -will be

attained, however, for any A with eigenvalues for which the corresponding t. include all the p

+ 1 values t with |T (t) | » 1. For such

i P

A, the x

q

for which [f(x-^)/f (xQ

) ] 2 = p are in the subspace spanned byjust those eigenvectors u^ belonging to the A

^ for which

XT ( A 3 | “ 1. The case p= 1 is exceptional, in that the maximum

jj^ is always attained in (3b) m (30), just because the two values

t - 1 where |T-^(t) j =* 1 are necessarily Z -images of eigenvalues.

For any fixed integer p

4 1 the above p-step process can be

iterated to yield a convergent procedure for finding A ^b. Namely,

for k “ 0, 1, •••

, one obtains as ihe unique point x in

VV for which f(x) is minimized. The optimum gradient method of

section £ is the case p

a 1, If p < < n we may expect that there

are eigenvalues A . of A close to the p + 1 points where
.

-

.
2 ?

| C XT ( A 3 | “ 1, and hence that the value of (3h) is almost attained for certain x . It is then to be expected that, for most x , 1

[f (xk+1)/f(x

k

)]

2 will be approximately

p

for all large k. Moreover,

if X

x << A n , then £ is near 1 and

p^

» 1 - 2p^( X A

n

) . Thus

the minimization in p dimensions may be expected asymptotically to

2 proceed p times as fast as the optimum gradient method (p 3 l) , when

2

p << n. When A /A-, “ n , as considered in section !?, the iterated

j.1 p-dimensional minimization in n-spaoe may be expected to converge like

the optimum gradient method in n/p dimensions.

The true asymptotic behavior of x^ - A is unknown. Does the

vector x^. - A ^b asymptotically lie in a certain (p + l) -dimensional

subspace, as is conjectured for p

= 1?

Because the above iterative process requires at each step solving

a linear system of p-th order, up to very recently the method was considered practical only for p << n. For p

* n, in particular, it

appeared that determining the minimizing o(

in (33) would involve

solving a linear system quite as difficult as Ax * b. Then, about

19^1, Stiefel [111;] , Lanczos [72], and Hestenes [!?0], working independently at first, all discovered that (33) can be minimized by

p repetitions of a beautifully simple algorithm. By taking p

- n

one actually finds A

_

^b in n steps, except for round-off errors. The

resulting conjugate gradient method is thus a typical finite iteration .

An extended exposition is in [5>.H] ? while brief expositions are given

by Rosser [100], Householder [£9], and Taussky and Todd [117],
' ,

: s

*

\

1
28

The conjugate gradient method is a nonlinear stationary iterative

process. The first approximant, is arbitrary one takes

p© “ r

Q

53 b •= Ax

q

. For ary k “ 0, assume the vectors x^., r

k

, p^ have

been determined,, Then r

k+ls ^k+1 are determined in order as

follows i

(35)

“k rW^k s Vi \

* “Vk 5 Vl - r

k

- « k

Apk

s \ “ -r

k+lAPk/PkAPk 5 Vi r

k*i

+ Vk •

Here r^ a b - Ax^ s and the significance of p^, °c^ 9 fly will appear

below. In the absence of round-off errors x

n

“ A ^b| if round-off

errors make x

R ^

A”^b 5 one has merely to carry the algorithm on for

k • n + 1, n + 2 a , until sufficient accuracy is attained.

The kernel of a number of methods of solving Ax - b for A > 0 is

the determination of a set of n directions (k • 0, •••

, n 1)

T which are conjugate (A-orthogonal) in the sense that p.Ap. 0 for

i / j . If the

£

p^j- are known,, then

(36) b E (Pkb/P^Apk

)pk

. k*o

A convenient method to get the is to apply the Gram-Schmidt process [99] of successively orthogonalizing some set of n linearly

29

independent vectors

|

. In Gaussian elimination (pivotal condensation) the v^. are taken to be the n coordinate unit vectors, as Fox,

Huskey, and Wilkinson [35] discovered, and the coefficients defining

the orthogonalization build up a triangular matrix.

In the method of conjugate gradients the v^ are the vectors

r

k

“ h - Ax^. (k s 0, , r - 1) • The beauty of this choice is that

rk+l turns out to be automatically conjugate to pQ

, p^, •••

, p^.

In picking P^+^_ it is therefore necessary only to alter r^.+ ^ in the

direction p^| the calculation of /3 and in (35) has this object.

The calculation of and is an iterative procedure for building

up the sum ( 36) .

Recalling that x actually minimizes (33), we see that in pracJr tice the residual r_, may become so small for some p < n that it is unnecessary even to complete the n steps theoretically required to compute

A ^b. This occurred for p « 90 in a successful calculation with n - 105

of a difficult stress problem reported by Hochstrasser [56] and Stiefel

• [111*] • Such a saving could hardly occur with Gaussian elimination,

since the unit vectors have no such intimate connection with the

system (l) as do the vectors r^.

The conjugate gradient method has a geometrical interpretation.

As in section one seeks the common center A of the ellipsoids

f(x) 88 constant. Let x

q

be prescribed. One gets x^ by performing

one optimum gradient step (section 5) in the unrestricted n-dimensional

space R . Recall the definition of ) . There is an (n - 1) -dimensional affine subspace R

n_-j_ passing through x^ and conjugate to ft"^(xo

)

.
.

,

'
30

The solution A lies in Rn=t^« One gets by taking an optimum

gradient step -within R^ The gradient p^ of f(x) at x^ within

R

n„i is the projection into R

n ^ of the gradient r.^ of f(x) at x^

within R^. The optimal point in the gradient direction p^ from x^ is

X

2

“ Similarly, one gets x^. +^ by taking one optimum gradient step from x^. within the (n - k)- dimensional affine subspace

Rn=k conjugate to iy^(x

o

) . Finally, R

q

is the solution

-1

point A b

.

This is a bare description of the method „ In [5b] Hestenes and

Stiefel give an amazing number of its properties, discuss its application to unsymmetrical systems (l), and so on. A few machine experiments [5l 9 56 s 5l 9 Hit] with the method suggest good stability with

/S.A. respect to round-off errors, but a theoretical study of the stability

would be desirable.

The conjugate gradient method can theoretically be applied to

solving a system Ax ® b, where A is a bounded positive-definite selfadjoint operator on a Hilbert space. One defines f(x) = (x. Ax)

- 2(x, b) + (b, A

=

1

b) j the inverse operator A~^ certainly exists.

The method will ordinarily no longer converge in n steps, but the asymptotic behavior of f(x ) can be discussed. Karush [ 70 ] shows that if

A «* oi I + K, where K is completely continuous and oc ^ 0, then f (x )

goes to 0 faster than the p-th term of any geometrical progression.

Hayes [U8] treats a general A with lower upper bounds

m, M (0 < m <M < 00 ), and proves that f(x^) “ [1 - (m/M) ]^f (x

q

)

.

More can be proved

;
'

'
31

v-1 Let S - (M + m)(M - m)”'

L

. The Birman inequality (3U) shows that

[f(x )/f(x )]

2 - l/T ( B) < 1. Hence , since 2T ( 8 ) - [£ + Vi> 2 - l]p

+ [S

(37 )

, one gets the estimate

2 f(x ) i/p

rs

r

_1

2p

s Vi2-i s +

In another paper the asymptotic nature of[f(x will be deP

scribed more precisely for a class of operators A with a continuous

spectrum.

8 . Errors and " condition .» One must say something about the

important but little-understood subject of errors. We may distinguish

between: I. errors committed in the course of solving the system by

a specific algorithm! H. errors inherent in the system Ax « b

.

Within I one is concerned with the truncation errors and the

round-off errors of an algorithm, a distinction explained in [12U]

The truncation error exists for infinite iterations, and may be

identified with x^ - A ^b! its behavior has been examined in the

above survey, under the assumption that there was no round-off error.

The study of round-off error itself is far more difficult, and there

seems to have been a complete discussion in connection with only one

method, elimination! see von Neumann and Goldstine [12U, h£]> and

also Mulholland [8E>]. For other studies see Bargmann, Montgomery,

and von Neumann [6], Dwyer [20], Satterthwaite [102], Thckennan

[121], and Turing [122].

/-W ^\_/V
-

'

*

-

!

'
32

Any approximate solution f of Ax ® b can be checked a posteriori

by forming the residual ? = b - A f . The magnitude of A""^b - I can

then be estimated by using some tool for examining errors under II.

Hence to bound the errors in a calculated c, it is unnecessary to

have an a priori knowledge of the accumulation of round-off error.

Such knowledge may be important, however, for planning purposes - for

example , in deciding in advance how many digital places to carry in

order that f be reasonably small.

The errors under H have attracted mere study. The practical

analyst, realizing that the elements of A and b are subject to uncertainty, wishes to know the corresponding uncertainty in A ^bj the

latter is the inherent error of Milne [81], The usual approach is the

approximate one of bounding the principal part of the error,

S (A

-1

b) ® A

"*( S A)A + A ^Sb; see Blumenthal [11] , Milne [81],

Moulton [8U], Ostrowski [90], Wittmeyer [126], and

Zurmuhl [131]. But others (e.g., Janet [61*], Ostrowski [89], Lonseth [??]) bound the entire error, finding a region S to which A ^b is

/w

rigorously confined, given f?

* b - A^ and other reasonably computable

quantities associated with the system. See also Woodbury [128],

Various persons (e.g., Jurgens [66], Todd [119], and Turing [122]}

have attempted to ascribe a condition to the system Ax b. In part,

the condition should measure the influence on A of small changes

in A and b| the larger the change in A "Hd for given changes in A and

b, the "worse" the condition. Although the condition depends on

both A and b, the measures hitherto proposed depend only on A. The
p.
33

belief is widespread that the condition of a system (l) has a decisive

influence on the convergence of an iterative solution and on the accuracy of a direct solution^ this cannot always be true. Even when it

is true for an iterative process, it may be possible actually to take

advantage of the poor condition of (l) in converting the slow process

into an accelerated method which converges rapidly. There is great

need for clarification of the group of ideas associated with "condition

With the concept of "ill-conditioned" systems Ax • b goes the idea

of "pre-conditioning" than. Gauss [111] and Jacobi [63] made early contrlbutions to this subject. That of Gauss is analyzed and extended in

[31].

4*/sJ A convenient means of pre-conditioning is to premultiply the

system with a matrix B, so that one has to solve

(38) BAx - Bb .

—1 T

The perfect choice of B would be A .A frequent choice is A , so

T

that (38) gets a symmetric matrix A A, very convenient for many processes, though "worse conditioned" in some senses (Taussky [UJ?] ) •

Is there in ary sense a "best choice" of B which is quickly obtainable?

Gauss’ elimination process may be written in the form (38), where

in the absence of round-off errors BA is a triangular matrix. In

some calculations the true solution A~^b comes from iteration of the

"back solution" - i.e., of getting x

k+1 by an approximate solution

of the triangular system BA(x - x^) - B(b - Ax^) . Where this occurs,

we may interpret the "forward solution" or triangularization as

merely a preconditioning of the system (l) into the form (38) •
; :

.. ,

*

'

. .

’

.

k
3h

9 • Influence of computing equipment. The usefulness of a process for solving (l) depends intimately on the properties of the

machine on which the calculation takes place,, as well as on the special

character of A and b. Hie past decade has seen revolutionary developments in computing equipment: analogue machinery, desk computers, IBM

equipment and automatically sequenced high-speed digital computers

As a result, computing methods are in no way settled down, and bibliographies are out of date before publication.

Analogue machinery can be very useful, but is not discussed herej

for references see Engineering Research Associates [22], Frame [36],

and Murry [86]

„

While direct methods for solving Ax ® b have been little mentioned here, they have been very successful since the time of Gauss or

earlier. Dwyer [20], Bodewig [12] and others conclude that a compact

arrangement of Gaussian elimination is commonly the best method for

the computing team of a desk machine, a data sheet, and a trained

human being - principally because the number of operations is minimized.

Elimination is very successful with IBM machines also, but its superiority

over other methods is less pronounced, because it is seldom expedient

with IBM equipment to use the compact arrangements which save so mary

operations. A bibliography on computing methods for IBM machinery is

given in [38a].

Let us now consider automatically programmed digital computers

like the SWAC. These are much faster than previous computers, but the

speed of arithmetic operations and of access to a small store of data
.

*

.

'
35

(high-speed memory) has been accelerated a great deal more than the

operations of input, output, and access to a large store of data.

The resulting change in the relative costs of different operations has

a profound effect on the choice of computing methods. O'ne soon learns

that a variety of processes have been tried and advocated for solving

(l)s certainly the optimal method depends on the problem, the machine,

the operator, and the coder. Moreover, small changes in these factors

may radically alter the optimal choice of method.

The following tentative ideas are based on a limited experience

with the ’SjiIAC, and practically no experience with other machines.

(The analysis is dominated by the relative shortage of memory cells in

the SWACs it is therefore less pertinent for machines with more storage

space, and for the 'SWAG after the expected addition of a magnetic

drum.) Assume n to be fairly large - say = 15. For simplicity we

again confine our attention to matrices A > 0. As indicated after

(38)

,

the forward solution in elimination amounts to building up a

new matrix BA, which must be stored somewhere. If BA is kept in the

high-speed memory, it occupies critically needed space. If it is

output and input as needed (say by rows, as Huskey [61] describes),

the programming is complicated and the solution is considerably

slowed. If A is a matrix which can be generated internally as needed

(for instance, the matrix of the Laplace difference operator) , it

requires little space, and BA becomes the principal item to be stored.

Where A cannot be generated internally, the storage problem gets

still worse, because the round-off errors can only be reduced by
.

'

>

.

-

'

'

’
36

using A to compute the residual r - b - Ax from time to time, so that

r both BA and A must be stored.

These considerations suggest that a solution process should preferably work with only one matrix, A itself, and should require relatively little other storage. Since the instructions have to be stored,

this suggests a process of simple structure, repeated as often as

necessary, A process seeming to require a minimum of storage is the

cyclic single-step procedure mentioned in section 3j besides A (if it

must be stored), one need store only one vector, x, at a time. This

method was picked by Reich [96] as best for a digital computer, and

is undoubtedly ideal when it converges fast enough. But we may expect

that the convergence is frequently too slow. If an acceleration of

the types discussed in section lj. is needed, the complication in programming may make another procedure preferable. Another method of

speeding up the cyclic single-step method is by appropriately overcorrecting (or undercorrecting ?) at each step, as discussed by

Frankel [38] and Young [129] for special systems (l) . It seans likely

that a careful determination of the optimal overcorrect!on will sometimes provide adequate convergence, but that it will often fail.

The ordinary relaxation (i,e,, piecewise linear) processes require

about the same storage as the cyclic single-step methods; it is not

clear whether they are essentially faster or not, A suggestion of

Motzkin and Schoenberg [83] for extreme overrelaxation is promising

but untried.

c
'

•
37

If the above methods fail, one can switch to the optimum gradient

method of section This also works with A, which must be stored or

generated, and further requires the storage of the two vectors x^. and

r^, (The storage of can be avoided if Xj^ - x^. is output at each

step, and cumulated later.) Again the method is probably commonly too

slow. It can be speeded tip either by the ft -optimum device of section

5», for ft < 1, or by Richardson's idea [98] of widely varying

*k - ^

a

k over the range of eigenvalues of A”"^,

If these tricks fail or require too complex a program, the grai dient methods of section 7 are available. Besides A, they require

the storage of three vectors x^, r^, p^» (As above, the outputting

of x^..,.^ “ saves storing x^.) Of these methods, there seems to be

no reason for not adopting the optimum gradient method, since for the

same storage its convergence is much the best. Programming is simple,

as only one routine is needed; all necessary variations in the

are provided automatically, A drawback is that, since the oc ^

and

other numbers vary so much in the calculation, It is difficult to

provide scaling factors in advance. Consequently one uses "floating

binaiy point" operations, requiring considerable memory space to

hold the instructions and multiplying the solution time by a factor

which varies on the SWAG from the order of £0 (when A is generated

internally) to about one (when A is input at each iterative step)

,

But the method has prcved able to cope with some very "badly conditioned" matrices, as reported by Hestenes, Hochstrasser, and

Wilson [E>1], and Hochstrasser [!?7]. It probably approaches the
'
38

ideal of a machine method which can be relied on to work automatically

without special analysis of the “particular system (l)

.

With any method the partitioning of A may greatly increase the

speed by enabling subsidiary matrix inversions to take place entirely

within the high-speed memory; see [21], One usually thinks of Gaussian

/'V'V elimination on the submatrix blocks. Would other methods on the

blocks be preferable?

39

BIBLIOGRAPHY

This is a selection of 133 titles out of over 5>00 in my card

file. It is obvious that many valuable references have had to be

omitted. Ordinarily, journal names are abbreviated according to

the list in Mathematical Reviews vol, 11 (I95fo) pp # 8£9-872 . Russian

is transliterated according tp the system loc. cit., p, 873 .

1. A,. A. Abramov, On a method of acceleration of iterative

processes (Russian), Doklady Akad. Nauk SSSR vol, Yh (19^0)

pp. 1051-10^2,

V' , 2, Shmuel Agmon, The relaxation method for linear inequalities , NAML Report £2-27, National Bureau of Standards, Los 'Angeles, 19f?l,

multilithbd typescript, 23 pp. To1 appear in Ann, of Math,

^3, A. C. Aitken, On Bernoulli’s numerical solution of algebraic

equations , Proc. Roy, Soc . Edinburgh, vol. U6 (1926) pp, 289“30£,

U. , Studies in practical mathematics. V. On the

iterative solution of a system of linear equations , Proc, Roy. Soc.

Edinburgh. Sect. A. vol. 63 ( 195> 0 ) pp„ £2-60 ?

D. G. Aronson: see [26],

£

.

T. Banachiewica, Methode de resolution numerique des

equations lineaires , du calcul des determinants et des inverses , et

de reduction des formes quadratiques Bull. Int. Acad. Folon. Sci.

Cl. Sci. Math. Nat. Ser. Sci, Math, (1938) pp. 393-HOlu Also in

Cracow observatory reprint 22 .
'

. <

,

*
Uo

6, V. Bargmann, D. Montgomery, and J, von Neumann, Solution of

linear systems of high order , Report prepared for the Bureau of

Ordnance (Contract NORD-9^96), 2!? October 19U6, multigraphed,. bound,

86 pp a

7, Commandant Benoit, Note sur une methode de resolution des

equations normales provenant de 1 * application de la methode des

moindres carres a un syst^me d" equations lineaires en nombre inf erieur

jt celui des inconnues . -— Application de la methode a la resolution

d* un systame defini d* equations lineaires (Procede du Commandant

Cholesky), Bull, Geodesique year 1921; (192I|.) no.. 2 pp 9 67"77o -

8 , M, D . Bingham, A new method for obtaining the inverse matrix, J, Amer. Statist, Assoc, vol, 36 (l9hl) pp. %30-H>3h9, M. ?. Birman, Some estimates for the method of steepest

descent (Russian), Uspehi Matem. Nauk (N,S,) vol, £ (19^0) no, 3

pp„ 1^2-1^ 0 Translated by C. D. Benster in NBS Report 2007,

National Bureau of Standards, Los Angeles, August, 195>2 .

10,

A. N. Black and R, V. Southwell, Relaxation methods applied

cJ~^S —"—* "

to engineering problems. II. Basic theory , with application to surveying and to electrical networks and an extension to gyros tatic systems .

Proc. Roy, Soc* London, Ser. A, vol, I6I4 (1938) pp, UU7-U67.

11* 0, Blumenthal, Uber die Genauigkeit der Wurzeln linearer

Gleichungen, Zeitschrift fur Math, und Physik vol, 62 (191U) pp, 3^9-362

12. E. Bodewig, Bericht uber die verschiedenen Methoden zur

Losung eines Sys terns line arer Gleiebongen mit reellen Koeffizienten.

I, II 3 HI * IV , V, Nederl, Akad, Wetensch,, Proc . vol, £0 (191*7)

pp. 930-9ljl, 110 Lr1116, 1285-1293, and vol a 5l (191*8)* pp. 53-6U,

and 211-219 = Nederl, Akad,; Wetensch. Indagationes Math. vol. 9 (191*7) -

pp. Ula-1^2, 518-530, 611-621, and vol n 10 (191*8). pp. 2U-35 and 82-90.

13. H. Boltz, Entwicklungsverfahren --- zur Ausgleichung geodatischer - ~

~

' " ,rr,r ““ " Netze nach der Methode der kleinsten Quadrate , Yeroff entlichen des

Preussischen Geodatischen Institute N. F. no. 90 (1923), Berlin.

li*. A* L. Cauchy, Methode gene rale pour la resolution des systemes

d 1 Equations simultanees 9 C. R. Acad. Sci» Paris vol, 25 (181*7)

pp. 536-538,

l5„ Lambertc Cesari, Sulla risoluzioni dei sistemi di equazioni

lineari per approssimazioni successive , Extract from Rassegna delle

poste, dei telegrafi e del teiefoni vol Q 1* (1937), 37 pp.

A. R. Collars see [39].

16. L, Collate. Graphische und numerische Verfahren, pp. 1-92 of P * omfrg »iiA ' ' v a ru.an,x.:i« ——— —— " • " ———»— Alwin Walther, Applied mathematics . Part I, FIAT Review of German

Science 1939*191*6, Office of Military Government, Wiesbaden, 191*8,

307 pp.

ft 17. — 9 uber die Konvergenzkriterien bei Iterationsverfahren

fur lineare Gleiehungs systeme

,

Math, Z, vol„ 5 3 (1950) pp. ll*9-l6l,

18. John H* Curtiss, S ampling methods applied to differential and

difference equations , pp . 87-IO9 of Proceedings , s eminar on scientific

computation November , 19h9 , International Business Machines Corp . , New

York, 1950, 109 pp
... .

,

.

.
19 . A . de la Garza An iterative method for solving systems of

linear equations , Report K-731, Carbide and Carbon Chem. Div., Union

Carbide and Carbon Corp., K-25* Plant,, Oak Ridge, Tenn., l£ Feb, 195>1,

mimeographed, l£ pp,

W. J, Duncans see [39].

20. Paul S. Dwyer, Linear computations , Wiley, New York, 19!?1,

3UU pp.

21. Eckert-Mauchly Div. of Remington-Rand Corp., Matrix algebra

programs for the UNIVAC

,

Philadelphia, 19^1, typescript, 11 pp.

22. Engineering Research Associates, High-speed computing devices , McGraw-Hill, New York, 19!? 0, U£l pp. (Supervisors: C. B, Tompkins and

J. H. Wakelinj editor? W. W, Stifler, Jr.)

23. 7*. N. Faddeeva, BHHHCiroTe'ji&HHe Me1 toah jmEfeiiHOH ajirebpu (computational methods of linear algebraj Russian), Moscow-Leningrad, 19^0,

2I4.O pp. (A translation is being prepared at the National Bureau of

Standards, Los Angeles.)

21;. Donald A. Flanders and George Shortley, Numerical determina-

'"Vv

’ ” ” tion of fundamental modes , J„ Appl p Phys . vol. 21 (19!?0) pp, 1326-1332 ?

2£. A. I, and G, E. Forsythe, Punched- card experiments with

accelerated gradient methods for linear equations , NBS Report 16U3,

National Bureau of Standards, Los Angeles, March 19!?2, multilithed

typescript, 29 pp. To appear in J, Research Nat, Bur. Standards,

26. George E, Forsythe, Theory of selected methods of finite

matrix inversion and decomposi tion. Lecture notes by D. G. Aronson and

K, E. Iverson^ INA Report 5>2-5>, National Bureau of Standards, Los

Angeles, 19!? 1, multilithed typescript, 93 pp.
' ' : r:Y\ v -r • ; . ’
27 , , Tentative «|.« rm n —~ai.na. ni " classification of methods and iat - - r r

.. . "" bibliography on solving systems of linear equations , INA Report 5>2-7,

National Bureau of Standards, Los Angeles, 19^1, multilithed typescript

78 pp, To appear in [87]. VON-* 28. George E, Forsythe and Richard A , Leibler, Matrix inversion by

a Monte Carlo method , Math. Tables and Other Aids to Computation vol, U

(19^0), pp. 127-129 and vol . 5 (19*1) p,

29 o G. E. Forsythe and T. S. Motzkin, Asymptotic properties of

the optimum gradient method , abstract, Bull, Amer. Math. Sac . vol, £7

(19£l) p. 183.

30. , Acceleration of the optimum

gradient method . Preliminary report , abstract, Bull, Amer. Math, Soc.

vol. £7 (19*1) pp, 30U-30£,

31. , An extension of Gauss » - — "f "~ 1 1 " transforma tion for improving the condition of systems of linear equations , Math t Tables and Other Aids to Computation vol, 6 (19^2),

« pp, 9-17.

32. L, Fox, A short account of relaxation' methods. Quart. J.

Mech. Appl, Math, vol, 1 (19U8) , pp. 2^3-280,

33. — , Practical solution of linear equations and inversion ' ri " ' “ T

"

of matrices , manuscript, about 60 pp. To appear in J . Research Nat.

Bur. Standards.

3h. Practical methods for the solution of linear equations and the inversion of matrices , J. Roy. Statist. Soc, Ser. B.

vol. 12 (1#0) * pp. 120-136.
.

.

... . •;

„

=

; ; .. -

, .

‘

. .

- i

;

, :

...... •:

• I

/ ' . . • .

-

.

’

- ;;.:A . -

:
uu

35. L. Fox, H„ D, Huskey, and J, H, Wilkinson, Notes on the solution of algebraic linear simultaneous equations . Quart. J. Mech. Appl.

Math, vol . 1 (19U8) pp. lIj.9-173 o

36. J. S. Frame, Machines for solving algebraic equations

,

Math.

Tables and Other Aids to Computation vol. 1 (191j5) pp. 337-353.

37 » [Stanley P. Frankel], Bibliography on computing machines,

/V>v

' —————— — Analysis Laboratory, Calif. Inst, of Technology, Pasadena, October 19149,

hectographed typescript, HI; pp. ( Frankel* s name not on copy.)

38. Stanley P. Frankel, Convergence rates of iterative treatments

of partial differential equations , Math. Tables and Other Aids to Computation vol. U (1950) pp. 65-?5.

ij8a« Ann D. Franklin and Eric V. Hankam, Bibliography on the

Business Machines Corp., New York, 1952, 50 pp. See also [62].

39. R, A, Frazer, W. J, Duncan, and A. R. Collar, Elementary

—V'V t matrices and some applications to dynamics and differential equations ,

Cambridge Univ. Press, 1938, Ul6 pp.

UO. Konrad Friedrich and Werner Jenne, Geometrisch-anschauliche s Auflosung linearer mit Nullkoeffizienten ausgestatteter Gleichungssys teme , Deutsche Akad. Wiss. Berlin. Veroff. Geodat. Inst. Potsdam,

no, 5 (l95l), 68 pp.

I4I. C. F. Gauss, Letter to Gerling, 26 December 1823, Werke ,

tionum erroribus minimis obnoxiae , 1826, Werke , Gottingen, vol. U

use of IBM machines in science , statistics, aid education . International

Gottingen, vol. 9 pp. 278-281. Translated by G. E. Forsythe,

Math. Tables and Other Aids to Computation vol. 5 (l95l)

pp f 255-258.

1|2 . , Supplementum theoriae combinationis observapp, 55-93.
.

’

'
1*3. M. K. Gavurin, Application of polynomials of best approximation to improving the convergence of iterative processes (Russian)

,

Uspehi Matem. Nauk (N.S.) vol. 5 (1950) no. 3 pp. 156-160. Translated by C. D, Benster in UBS- Report 2007, National Bureau of

Standards, Los Angeles, August 1952.

U*. Hilda Geiringer , On the solution of systems of linear equations by certain iterati on methods , pp. 365-393 of Reissner Anniversary

Volume, Contributions to applied mechanics . Edwards, Ann Arbor, Mich.,

191*9.

s see also [ 123] .

1*5. Herman H. Goldstine and John von Neumann, Numerical inverting of matrices of high order, II , Proc. Amer. Math. Soc. vol, 2

(1951) pp. 188-202. Part I is [121*3.

Eric V. Hankaras see [38a],

1*6. D. R. Hartree, Experimental arithmetic , Eureka vol. 10 (191*8)

pp. 13-18,

1*7 • Harvard Computation Laboratory, A manual of operation for the

automatic sequence controlled calculator. Harvard Univ, Press, 191*6,

561 pp,

1*8. R, M. Hayes, Iterative methods of solving linear problems on

Hilbert space , Ph.D. Thesis, Univ. of California, Los Angeles, 1952,

multilithed typescript, 62 pp.

1*9. A, Hertwig, Die Losung linearer Gleichungen durch unendliche

Reihen und ihre Anwendungen auf die Berechnung hochgradig statisch

unbestimmter Systeme , pp. 37“59 of Festschrift fur H. Muller-Breslau,

Leipzig, 1912.
•

'

.

'

'
U6

£G„ M. Ha Hestenes, Iterative methods for solving linear equations , NAML Report E>2-9, National Bureau of Standards , Los Angeles, 19!?1,

multilithed typescript, 19 pp.

5>1* M, R. Hestenes, Urs Hochstrasser, and L. S. Wilson, Some

numerical examples on solving systems of linear equations by the conjugate gradient method for non-symmetric systems of equations , National

Bureau of Standards, Los Angeles, manuscript.,

f>2» Magnus R. Hestenes and William Karush, A method of gradients

for the calculation of the characteristic roots and vectors of a real

symmetric matrix , J. Research Nat. Bur. Standards vol* Ii7 (l95>l)

pp. Hj-61.

f>3* Magnus R. Hestenes and Marvin L. Stein, The solution of linear

equations by minimization , NAML Report S>2-lj£, National Bureau of Standards, Los Angeles, 12 December 195>1, multilithed typescript, 35> pp*

5>lu Magnus R„ Hestenes and Eduard Stiefel, Method of conjugate

*v-\^ gradients for solving linear systems, NBS Report 16!? 9 s National Bureau

a of Standards, Los Angeles, 10 March 195>2, multilithed typescript,

101 pp . To appear in J. Research Nat. Bur. of Standards,

3>E>. Thomas J. Higgins, A survey of the approximate solution of

n/'v *— —

— ' —

two-dimens ional physical problems by variational methods and finite

difference procedures , pp e 169-198 of L. E. Grinter (editor),

Numerical methods of analysis in engineering , MacMillan, New York,

19U9, 207 PP*
.

.

.
56. £ Ubrs 3 Hochstrasser, Die Berechnung der Verschiebungen in

einer parallelograimnformigen Scheibe , Bericht S-17, Techn. Abteilung,

Eidgenossisches Flugzeugwerk, Emmen, Switzerland, rani tigraphed typescript, 21 pp.

57. Urs Hochstrasser, The solution of systems of linear equations

& the conjugate gradient method for use on IBM equipment , National

Bureai of Standards, Los Angeles, manuscript*

___

s see also [5l]

,

58. Harold Hotelling, Some new methods in matrix calculation ,

Ann. Math, Statistics vol, lU (19U3) pp. 1-33 and I4J4D—UUl

59 0 A. S. Householder, Notes on numerical methods , Oak Ridge

National Laboratory, Tenn., 19b9-1952, multilithed typescript,

c. hOO pages so far.

60*

Some numerical methods for solving systems

of linear equations , Amer. Math, Monthly vol. 57 (1950) , pp. Ii53-U59.

6l o H„ D. Huskey, Characteristics of the Institute for Numerical

Analysis computer , Math. Tables and Other Aids to Computation vol. U

(1950), pp* 103-108.

s see also [35].

62 r International Business Machines Corp., Bibliography on the

use of IBM machines in science , statistics , and education . New York,

1950 „ Apparently superseded by [38a].

K« E e Iverson? see [26] .
*

'

.
U8

63. C. G. J. Jacobi, Ueber eine neue Auflo sungs art der bei der

Metho de der kleinsten Quadrate vorkommenden linearen Gleichungen,

Astr. Nachr. vol. 22 (l8Ij5>) no. 323 pp. 297-306. (Also in Jacobi’s

Werke vol. 3- p. U6??)

6k* M. Janet, Sur les syst^mes d' Equations aux derivees partielles,

C. R. Acad. Sci, Paris vol. 170 (1920 ) pp , 1101-1103 and J. Math.

Pures Appl • Ser. 8, vol, 3 (1920), pp. 63-l3l.

Werner Jen ne 1 see [I4O].

63 • Henry Jensen, An attempt at a systematic classification of

seme methods for the solution of normal equations , Meddelse No. 18,

Geodaetisk Institut, Copenhagen, 19bb„ b$ pp.

66. Enno Jurgens, Zur Auflosung linearer Gleichungssysteme und

numerischen Berechnung von Determinanten , Palm (printer), Aachen, 1886,

20 pp,

«• •» 67* S, Kaczmarz, Angenaherte Auflosung von Sjystemen linearer

Gleichungen , Bull. Int. Acad. Polon. Sci. Cl, Sci. Math. Nat. Ser, A.

Sci. Math. (1937) pp. 333-337.

68. L, V. Kantorovic, Functional analysis and applied mathematics

(Russian), Uspehi Matem. Nauk (N.S.) vol. 3 (19U8) no . 6 pp, 89-183.

(A translation is being issued by the National Bureau of Standards,

Los Angeles .)

, v

69. L, V. Kantorovic and V. I y Krylov, IIpHdjrHSceHHHe ”'“rt weTOflH bhc- " 1 1

urero aHaJiH3a (Approximation methods in higher analysis; Russian),

3rd edit., Moscow-Leningrad, 1930, 693 pp.

U9

JO. W* Karush, Convergence of a method of solving linear problems , NAML Report 52-19 , National Bureau of Standards, Los Angeles,

21 August 1951, multilithed -{Typescript, 20 pp. To appear in Amer. J.

Math. — : see also [52] .

70a. M. A. Krasnosel^ skii' and S. G. Kre^n, The iterative process

with minimal residuals (Russian), Mat. Sbornik (N.S.) vol. 31 (73)

(1952) pp. 315-33H.

V. I. Krylovs see [69].

71. A. G. Kuros, A. I. Maricusevic, and P. K. Ra^evskx?. MaxeMaTHKa

b_ GGCP 3 a TpwflpaTb jie t 1917-19U7 (Mathematics in the USSR in the tiiirty

years 19l7-19ii7$ Russian), Moscow-Leningrad, 19U8 , IOI4I4. pp. (Pp. 75 9-857

on numerical and graphical methods contain a survey and comprehensive

bibliography of Tfcussian work.)

72. Cornelius Lanczos, Solution of —systems —— of linear equations by ' mmJLm minimized iterations , J. Research Nat. Bur. Standards vol. U9 (1952)

pp. 33-53.

73. , Chebyshev polynomials in the solution of

large scale linear systems . National Bureau of Standards, Los Angeles,

manuscript

.

7h. H. Liebmann, Die angenaherte Ermittlung harmonischer

/VV

Funktionen und konformer Abbildung, Sitzungsberichte Bayer. Akad.

Wiss. Math .-Nat. Abt. Math.-Phys. Kl. vol. U7 (1918) pp. 385-I|l6.

75, A. T. Lonseth, The propagation of error in linear problems,

Trans. Amer. Math. Soc. vol. 62 (19U7) pp. 193-212.
1 ' '. , ' - v. • ' • '> ;
So

76. Samuel Lubkin,‘ A method of summing infinite' series , J. Research ”” ————— — r " “ ’ r Nat. Bur. Standards vol. U8 (195>2) pp. 228-2SU.

77. L. A. Lyusternik, Remarks on the numerical solution of boundary problems for Laplace* s equation and the calculation of characteristic values by the method of nets (Russian), Trudy Mat. Inst. Steklov

vol* 20 (19U7) pp. U9-6U.

78. C. C. MacDuffee, The theory of matrices , Chelsea, New York,

19U6, 110 pp ?

79. W, Markoff, Uber Polynome , die in einem gegebenen Intervalle

moglichst wenig von Null abweichen, Math. Ann. vol. 77 (1916)

pp. 213-2^8. (Translation and condensation by J . Grossman of Russian

article published In 1892 .)

A. I. Markusevics see [71] .

80 „ S . N. Mergelyan, Uniform approximation of functions in the

complex plane (Russian), Uspehi Matem. Nauk (N.S.) vol. 7 (195>2) no. 2,«

pp, 31-122. Translated in Amer. Math. Soc. Translation T-113>

.

81, William Edm.und Milne, Numerical calculus, Princeton Uhiv,

Press, 19U9, 393 pp.

D. Montgomery? see [6],

82. Theodore S. Motzkin, Bibliography on linear inequalities,

linear programming, game strategy, economic behavior, and statistical

decision functions, in preparation for probable issue by the National

Bureau of Standards, Los Angeles.

see a^- so [29 , 30, 31> 83].
*

.

'
5l

83. T. S. Motzkin and I» J. Schoenberg, to. the relaxation method

for linear inequalities , MBS Report 1881, National Bureau of Standards,

Los Angeles, August 1952, multilithed typescript, 21 pp. To appear in

Proc. Araer. Math. Soc.

8U. F. R. Moulton, On the solutions of linear equations having

small determinants , Amer. Math. Monthly vol. 20 (1913) pp. 2U2—2U9

85 . H , P • Mulholland , On the distribution of a convex even function of several independent rounding-off errors , Proc. Amer. Math.

Soc. vol. 3 (1952) pp. 310-321.

86. Francis J . Murray, The theory of mathematical machines , revised edition, King»s Crown Press, New York, 19U8, 139 pp,

87 . National Bureau of Standards, Proceedings of a symposium on w

simultaneous linear equations and the determination of eigenvalues,

August 195l 0 In press.

88. P. A c Nekrasov, Deteimination of the unknowns by the method

of least squares for very many unknowns (Russian), Mat. Sbornik

vol, 12 (188U) pp. 189-201*.

89. Alexandre Ostrowski, Sur la determination des homes

inferieures pour une class e des determinants , Bull. Sci . Math. Ser. 2 r

vol. 61 (1937) pp. 19-32.

90 . f Sur la variation de la matinee inverse

d'une matrice donnee , C. R„ Acad. Sci. Paris vol. 231 (1950)

pp. 1019-1021

£2

91. , Determinants with dominant principal W

diagonal and the absolute convergence of linear iteration processes

(text in German) , NBS Report I727, National Bureau of Standards,

Washington, 1 June 19!?2, multilithed typescript, U9 pp.

92^. 9 to the convergence of cyclic linear iterations for symmetric and nearly symmetric matric es . U, NBS Report 17!?9 j National Bureau of Standards, Los Angeles, 26 June 19!?2, multilithed

typescript, £ pp.

93 • ______ ___________ s to. '^ie linear iteration procedures for

symmetric matrices , NBS Report I8I4I4., National Bureau of Standards,

Los Angeles, I; August 195>2, multilithed typescript, 29 pp.

9iu Sam Perlis, Theory of matrices , Addison-Wesley, Cambridge,

Mass,, 19J>2, 237 pp,

9!>. P. Pizzetti, Sulla eompensazione delle osservazioni secondo

il metodo dei minimi quadrati , nota I, II, • Rendiconti della Reale

Accademia Nazionale dei Lincei, Rome, Ser. b9 vcl, 3 (1887) pp. 230-23!?

and 288-293.

Hilda Pdllaczek-Geiringers see [UU] and [123].

A—

Pj K. Rasevskiis see [71].

96. Edgar Reich, The solution of linear algebraic equations by

successive approximations , Memorandum M-969 , Servomechanisms Laboratory, Mass. Inst, of Technology, Cambridge, 5> August 19U8, hectographed, 36 pp.

97 , , On the convergence of the classical iterative

method of solving linear simultaneous equations , Ann, Math, Statistics

vol • 20 (19U9) pp p UU8-l£l.
_
S3

98. L. F. Richardson, The approximate arithmetical solution by

finite differences of physical problems involving differential equations , with an application to the stresses in a masonry dam, Philos

.

Trans. Roy. Soc. London Ser. A. vol. 210 (1910) pp. 30?-35>7.

99. Maria Sofia Roma, Sulla risoluzione numerica dei sistemi di

equazioni algebriche linear! eol rnetodo della ortogonalizzazione , Ricerca Sci. vol. 20 (19!?0) pp. 1288-1290 * Consiglio Naz. Ricerche.

Pubbl. 1st. Appl. Calcolo no. 283.

R. L. Rosenbergs see [113].

100^. J . Barkley Rosser, Rapidly converging iterative methods for

solving linear equations „ To appear in [ 87 ].

101. Paul A. Samuelson, A convergent iterative process, J. Math.

Physics vol. 2\± (l9li>) pp. 131-13U.

102^, F. E. Satterthwaite, Error control in matrix calculation ,

Ann. Math. Statistics vol„ 15> (19UU) PP . 373-387.

I. J. Schoenbergs see [ 83 ].

103. E. Schroder, Uber unendlich viele Algorithmen zur Auflosung

der Gleichungen, Math. Ann„ vol. 2 (I87O) pp. 3X?v36!?„

10U. G. Schulz, Iterative Berechnung der rezip roken Matrix ,

Z. Angew. Math. Mech. vol. 13 (1933) pp. $1~69 •

10£. I. Schur, Uber Potenzreihen , die im Innern des Einheitskreises

beschrankt sind, J. Reine Angew Math. vol. llj.7 (1917) PP» 20^-232.

106. [Hans Schwerdtfeger] , Bibliography on iteration , [reproduced

at National Bureau of Standards, Los Angeles, 19!? 1] , multilithed

typescript, 13 pp. (Bracketed material not on copy.)
.
107. Ludwig Seidel , Ueber ein Verfahren , die Gleichungen, auf

welehe die Methode der kleinsten Quadrate fuhrt , sowie lineare

o© ab — Gleichungen uberhaupt , durch successive Annaherung aufzulosen , Abh.

Bayer. Akad. Wiss . Math .-Nat. Abt. vol. 11 (I87U) no . 3 .pp . 81-108.

108. Daniel Shanks, An analogy between transients and mathematical

sequences and some nonlinear sequence-to-sequence transforms suggested

by it. Part I, NOLM 999U, Naval Ordnance Lab'., Silver Spring 19, Md.,

26 July 19h9, multilithed typescript, U2 pp.

Short!ey5 George: see^^].

109. R. V. ' Southwell, Relaxation methods in engineering science, a treatise on approximate computation , Oxford Univ. Press, 19i|0, 5

2!?2 pp,

110. , Relaxation methods in theoretical physics ,

Oxford Univ, Press, 19U6, 21$ pp.

: see also [10].

111. J. L. Steam, Iterative solutions of normal equations , Bull.

Geodesique year 195uL (195>1) PP- 331-339.

112. Marvin L. Stein, Gradient methods in the solution of systems

of linear equations, J. Research Nat. Bur. Standards vol. U8 (19!?2) ~

pp. I4.O7—Ijl3 .

(/)

0 )

% see also [f>3]

.

113.

P. Stein and R. L. Rosenberg, On the solution of linear

simultaneous equations by iteration , J. London Math. Soc. vol. 23

(19U8) pp. 111-118.

11U. E . Stiefel, Uber einige Methoden der Relaxationsrechnung ,

Z. Angew. Math. Physik vol. 3 (19E>2) pp. 1-33.

s see also [9U]

.
* * ' •
W. W. Stifler, Jr.; see [22],

/Vv» ll£. 0. Taussky, Note on the condition of matrices , Math. Tables

and Other Aids to Computation vol . I; (1950) pp. 111-112.

116. Olga Taussky, Bibliography on bounds for characteristic roots

of finite matrices , MBS Report 1162, National Bureau of Standards,

Washington, September 195>1, multilithed typescript, 10 pp.

117. Olga Taussky and John Todd, Systems of equations , matrices

and determinants, Math. Mag. vol, 26 (19^2) pp. 9-20 and 00-00.

118 . G . Temple, The general theory of relaxation methods applied

to linear systems , Proc. Roy, Soc. London Ser. A. vol. 169 (1939)

pp. U76-£00 .

119. John Todd, The condition of a certain matrix , Proc. Cambridge

Philos. Soc. vol. U6 (19U9) pp. 116-118.

• see also [ 117 ]

.

120. C. Tompkins, Proj ection methods in calculation of some linear

problems , [19U9?], multilithed typescript, £2 pp., part of Engineering

Research Associates, Logistics papers , issue no . IV, Appendix I to

Bimonthly Progress Report No. 17, Contract N6onr-2l;0, Task 0'rder I,

Office of Naval Research, Project NR 0l;7 010. —— —i see also [22],

121. L. B. Tuckerman, On the mathematic ally significant figures

in the solution of simultaneous linear equations , Ann. Math. Statistics vol. 12 (l9lil) pp. 307-316.

122. A. M. Turing, Rounding-off errors in matrix processes,

Quart. J. Mech. Appl. Math. vol

? 1 (19U8) pp. 287-308.
.

.

•

fvjr =;• * .

.

.

. .

*

.
£6

123# R. von Mises and Hilda Pollaeask-Geiringer, Praktische

Verfahren der Gleichungsauflosung , Z. Angew. Math. Mech. vol. 9 (1929)

pp. £8-77 and l£2~l6h.

J. von Neumann: see [6] and [!}£].

/V* /w

12U. John von Neumann and H # H. Goldstine, Numerical inverting of

matrices of high order. Bull. Amer. Math. Soc. vol, £3 (19U7)

pp. 1021-1099. Part II is [lj£].

J. H. Wakelin: see [22].

12£ . Wolfgang Wasow, A note on the inversion of matrices by random

•walks . Math. Tables and Other Aids to Computation vol. 6 (19£2)

pp. 78-81.

J. H, Wilkinson: see [3£] .

L. S. Wilson: see [£l] .

•• 126. Helmut Wittmeyer, Einfluss der Anderung einer Matrix auf die

Losung des zugehorigen Gleichungssystems , sowie auf die characteristischen Zahlen und die Eigenvektoren, Z. Angew. Math. Mech, vol. 16

(1936), pp. 287-300.

•® •• 127. , Uber die Losung von linearen Gleichung

s

-

systemen durch Iteration , Z. Angew. Math. Mech. vol. 16 (1936)

pp. 301-310.

128. Max A. Woodbury, Inverting modified matrices . Memorandum

Report h2. Statistical Research Group, Princeton, ll|. June 19£0,

hectographed, H pp.

£7

129. David M. Young, Jr., Iterative methods for solving partial

difference equations of elliptic type , Ph.D. Thesis, Harvard Uhiv.

Math. Dept., 1$^>0, blueprint, c. 100 pp,

130 . S . I . Zuhovicki'i, An algorithm for the solution of the

V y Cebysev approximation problem in the case of a finite system of incompatable linear equations (Russian), Doklady Akad„ Nauk SSSR vol. 79

(19^1) pp. £6l-£6U.

131* Rudolf Zurmuhl, Matrizen . ELne Darstellung fur Ingenieure ,

Springer, Berlin-GSttingen-Heidelberg, 19^0, I4.27 pp.

National Bureau of Standards, Los Angeles

November 18, 19^2

THE NATIONAL BUREAU OF STANDARDS

Functions and Activities

The functions of the National Bureau of Standards are set forth in the Act

of Congress, March 3, 1901, as amended by Congress in Public Law 619, 1950.

These include the development and maintenance of the national standards of

measurement and the provision of means and methods for making measurements

consistent with these standards; the determination of physical constants and

properties of materials; the development of methods and instruments for testing

materials, devices, and structures; advisory services to Government Agencies on

scientific and technical problems; invention and development of devices to serve

special needs of the Government; and the development of standard practices,

codes, and specifications. The work includes basic and applied research, development, engineering, instrumentation, testing, evaluation, calibration services

and various consultation and information services. A major portion of the

Bureau’s work is performed for other Government Agencies, particularly the

Department of Defense and the Atomic Energy Commission. The scope of activities is suggested by the listing of divisions and sections on the inside of the

front cover. Reports and Publications

The results of the Bureau’s work take the form of either actual equipment and

devices or published papers and reports. Reports are issued to the sponsoring

agency of a particular project or program. Published papers appear either in the

Bureau’s own series of publications or in the journals of professional and scientific

societies. The Bureau itself publishes three monthly periodicals, available from

the Government Printing Office: The Journal of Research, which presents complete papers reporting technical investigations; the Technical News Bulletin, which

presents summary and preliminary reports on work in progress; and Basic Radio

Propagation Predictions, which provides data for determining the best frequencies

to use for radio communications throughout the world. There are also five series

of nonperiodical publications: The Applied Mathematics Series, Circulars, Handbooks, Building Materials and Structures Reports, and Miscellaneous Publications.

Information on the Bureau’s publications can be found in NBS Circular 460,

Publications of the National Bureau of Standards ($1.00). Information on cali- bration services and fees can be found in NBS Circular 483, Testing by the National

Bureau of Standards (25 cents). Both are available from the Government Printing Office. Inquiries regarding the Bureau’s reports and publications should be

addressed to the Office of Scientific Publications, National Bureau of Standards,

Washington 25, D. C.