pragma License (Unrestricted);

AI95-00296

!summary

The vector and matrix operations in ISO/IEC 13813 plus related operations are added to Ada.Numerics in Annex G.

!problem

A number of secondary standards for Ada 83 were produced for the numerics area. Most of the functionality of these standards was incorporated into Ada 95 (some in the core language and some in the Numerics Annex) but two packages from ISO/IEC 13813 were not. These are generic packages for the manipulation of real and complex vectors and matrices.

The UK was asked to review the situation and to make a recommendation; they recommended that if Ada were amended then consideration should be given to including the packages within the Numerics annex.

The packages can be implemented entirely in Ada and thus present little burden to implementors. Providing secondary standards has not proved satisfactory because they are not sufficiently visible to the user community as a whole.

!proposal

It is proposed that two generic packages be added to the numerics annex. They are Ada.Numerics.Generic_Real_Arrays and Ada.Numerics.Generic_Complex_Arrays. They are included as a new subclause G.3 in order to avoid excessive renumbering of other clauses.

In addition to the main facilities of 13813, these packages also include subprograms for the solution of linear equations, matrix inversion, determinants, and the determination of the eigenvalues and eigenvectors of real symmetric matrices and Hermitian matrices.

The library package Numerics.Real_Arrays is declared pure and defines the same types and subprograms as Numerics.Generic_Real_Arrays, except that the predefined type Float is systematically substituted for Real'Base throughout. Nongeneric equivalents for each of the other predefined floating point types are defined similarly, with the names Numerics.Short_Real_Arrays, Numerics.Long_Real_Arrays, etc.

Two types are defined and exported by Ada.Numerics.Generic_Real_Arrays. The composite type Real_Vector is provided to represent a vector with components of type Real; it is defined as an unconstrained, one-dimensional array with an index of type Integer. The composite type Real_Matrix is provided to represent a matrix with components of type Real; it is defined as an unconstrained, two-dimensional array with indices of type Integer.

The effect of the various functions is as described below. In most cases the functions are described in terms of corresponding scalar operations of the type Real; any exception raised by those operations is propagated by the array operation. Moreover the accuracy of the result for each individual component is as defined for the scalar operation unless stated otherwise.

In the case of those operations which are defined to involve an inner product, Constraint_Error may be raised if an intermediate result is outside the range of Real'Base even though the mathematical final result would not be.

Implementation Requirements

Accuracy requirements for the subprograms Solve, Inverse, Determinant, Eigenvalues and Eigensystem are implementation defined.

For operations not involving an inner product, the accuracy requirements are those of the corresponding operations of the type Real in both the strict mode and the relaxed mode (see G.2).

For operations involving an inner product, no requirements are specified in the relaxed mode. In the strict mode the modulus of the absolute error of the inner product X*Y shall not exceed g*abs(X)*abs(Y) where g is defined as

   g = X'Length * Machine_Radix**(1-Machine_Mantissa)

Implementations shall document any techniques used to reduce cancellation errors such as extended precision arithmetic.

AARM Note

The above accuracy requirement is met by the canonical implementation of the inner product by multiplication and addition using the corresponding operations of type Real'Base and performing the cumulative addition using ascending indices. Note however, that some hardware provides special operations for the computation of the inner product and although these may be fast they may not meet the accuracy requirement specified. See Accuracy and Stability of Numerical Algorithms By N J Higham (ISBN 0-89871-355-2), Section 3.1.

Implementation Permissions

The nongeneric equivalent packages may, but need not, be actual instantiations of the generic package for the appropriate predefined type.

Implementation Advice

Implementations should implement the Solve and Inverse functions using established techniques such as LU decomposition with row interchanges followed by back and forward substitution. Implementations are recommended to refine the result by performing an iteration on the residuals; if this is done then it should be documented.

It is not the intention that any special provision should be made to determine whether a matrix is ill-conditioned or not. The naturally occurring overflow (including division by zero) which will result from executing these functions with an ill-conditioned matrix and thus raise Constraint_Error is sufficient.

The test that a matrix is symmetric may be performed by using the equality operator to compare the relevant components.

!discussion

Section G.1.1 of the Rationale for Ada 95 says

"A decision was made to abbreviate the Ada 95 packages by omitting the vector and matrix types and operations. One reason was that such types and operations were largely self-evident, so that little real help would be provided by defining them in the language. Another reason was that a future version of Ada might add enhancements for array manipulation and so it would be inappropriate to lock in such operations prematurely."

It is now clear that such enhancements will not be added so the second justification for omitting the facilities of 13813 disappears. In order to overcome the objection that everything is self-evident we have taken the approach that we should add some basic facilities that are commonly required, not completely trivial to implement but on the other hand are mathematically well understood.

The overall goal is thus twofold

* to provide commonly required facilities for the user who is not a numerical professional,

* to provide a baseline of types and operations that form a firm foundation for binding to more general facilities such as the well-known BLAS (Basic Linear Algebra Subprograms, see www.netlib.org/blas).

The packages closely follow those in 13813. However, the discussion has been considerably simplified by assuming the properties of the corresponding scalar operations such as those in Numerics.Complex_Types. This removes a lot of explicit mention of raising exceptions and accuracy. Also remarks that these are standard mathematical operations have been deleted when the meaning is given by other words in the description.

The component by component operations of * / and ** on vectors have been deleted on the grounds that they are not useful. (They might be useful for manipulating arrays in general but we are concerned with arrays used as vectors for linear algebra.)

Operations for vector products were considered but not added. This is because, as usually defined, they only apply in three-dimensional space.

It is hoped that there is not too much confusion between component when applied to the parts of a complex number and component when applied to a part of an array. 13813 uses element for the latter but the proper Ada term for an element of an array or field of a record is of course component.

Observe that the index range of the various arrays is Integer (rather than Natural or Positive). This permits negative indices and in particular arrays with symmetric index ranges about zero.

The function Identity_Matrix of 13813 has been changed to Unit_Matrix on the grounds that the prime concern is with linear algebra and not group theory.

The accuracy of most simple operations follows from the underlying operations on scalar components. In the case of inner product there is the potential for special operations to improve the speed and/or accuracy. We have specified reasonable requirements in the exact mode, which are met by the canonical implementation using a loop statement. This is because on some hardware, built-in instructions which are fast actually lose accuracy. Note that the Fortran language standard recognizes the existence of inner product as a special case but imposes no accuracy requirements at all.

These requirements for inner product are based on the analysis in Chapter 3 of Accuracy and Stability of Numerical Algorithms by N J Higham. A factor of nearly 2 has been added as a precaution. In the complex case the factor is increased by a factor of sqrt(2) in the mixed complex and real case to take account of the hypotenuse effect and by a further factor of 2 in the pure complex case since the number of terms is effectively doubled.

Functions have been added to Numerics.Generic_Real_Arrays for matrix inversion and related operations. On the grounds that it would be helpful to have simple algorithms available to users rather than none at all, no accuracy requirements are specified. Instead the accuracy is stated to be implementation-defined which means that the implementation must document it.

The names chosen are Solve and Inverse. The former is overloaded, one version solves a single set of linear equations; the second solves multiple sets. Note that Inverse is not strictly necessary because the effect can be obtained by using Solve and a Unit_Matrix thus

  I := Unit_Matrix(A'Length); B := Solve(A, I); -- same as B := Inverse (A);

A common technique when solving sets of linear equations is to refine the result by an iteration on the residuals. Thus to solve the set of equations Ax = y we first perform

  X := Solve(A, Y);

we then compute the error D by multiplying back thus

  D := Y - A*X;

and then compute the refinement to X

  DX := Solve(A, D);

and then correct X thus

  X := X + DX;

Implementations are recommended to do this automatically; it requires little extra computation if LU decomposition is used internally. If they do such refinement then it should be documented.

A function for computing the determinant has been added since it is of some help in deciding whether an array is ill-conditioned and therefore the results of inversion might be suspect. Determinants are also useful in some statistical calculations. The evaluation of determinants is very liable to overflow and many such routines return a scaling power of 10 in order to keep the basic result within range. For simplicity, it was decided not to do this since it is less likely with matrices of low order; the user can of course scale the components of the matrix if necessary.

Similar functions have also been added for complex arrays. However, it was not deemed necessary to provide for mixed real and complex operands for Solve.

In addition, subprograms have been added for the computation of eigenvalues and vectors of real symmetric matrices and Hermitian matrices. The subprograms are Eigenvalues and Eigensystem. There is no separate subprogram for eigenvectors since it is unlikely that these would be required without the eigenvalues.

The most common application is with real symmetric matrices and these have real eigenvalues and vectors. Applications include: moments of inertia, elasticity, confidence regions and so on.

A slight problem arises in deciding who should check that a matrix is symmetric, the user or the package. Computations of, for example, a covariance matrix will result in a matrix that ought to be exactly symmetric but small errors might be introduced which mean that it is not exactly symmetric. The onus could be placed on the user to ensure that the matrix is exactly symmetric or alternatively some tolerance could be passed to the eigen subprograms so that they can perform a reasonable check. Passing a tolerance level adds an irritating parameter, raises the issue of how the tolerance should be expressed and gives the user the problem of what tolerance to request. Moreover, the algorithms still have to decide which actual values should be used for those pairs of components that are not exactly equal.

Accordingly, we have placed the onus on the user to ensure that the matrix is exactly symmetric. This could be done for example by taking the mean of the matrix and its transpose. The test in the subprograms can then test for exact equality which is guaranteed to deliver the correct answer. If the test fails then Argument_Error is raised. This exception, declared in Ada.Numerics, is the normal exception for numeric operations when the argument is out of range. Note that errors such as when bounds of arrays do not match raise Constraint_Error by analogy with built-in array operations.

A third approach considered was for the user to supply a parameter indicating which half of the matrix should be used to define it (the upper or lower triangle). This avoids the need for any testing but it was considered bad practice to be able to pass junk in the other half of the matrix.

The complex equivalent of real symmetric matrices are Hermitian matrices. Hermitian matrices are such that their transpose (that is with rows and columns interchanged) equals their complex conjugate (that is with the sign of imaginary parts reversed). Hermitian matrices also have real eigenvalues. Applications include Quantum Mechanics.

Again we have placed the onus on the user to ensure that the matrix is Hermitian. The check in the package can then use strict equality for the real parts and negation followed by equality for the imaginary parts.

We considered providing subprograms for the determination of eigenvalues and eigenvectors of general real and complex matrices. Such matrices can have complex eigenvalues and therefore provision for these would have to be in the complex package. However, there are mathematical difficulties with these general cases which are in strong contrast to the real symmetric and Hermitian matrices. Thus, Numerical Recipes by Press, Flannery, Teukolsky and Vetterling says regarding the real case:

  "The algorithms for symmetric matrices ... are highly satisfactory in
   practice. By contrast, it is impossible to design equally satisfactory
   algorithms for the nonsymmetric case. There are two reasons for this.
   First, the eigenvalues of a nonsymmetric matrix can be very sensitive
   to small changes in the matrix elements. Second, the matrix itself can
   be defective so that there is no complete set of eigenvectors. We
   emphasize that these difficulties are intrinsic properties of certain
   nonsymmetric matrices, and no numerical procedure can cure them."

Similar remarks apply to complex matrices where Hermitian matrices are well-behaved but non-Hermitian matrices can be troublesome.

In view of these computational difficulties and the fact that requiring the eigensystem of general matrices is uncommon, we decided not to provide such facilities.

Consideration was also given to the inclusion of explicit facilities for LU decomposition (as provided for example in the BLAS). LU decomposition is a common first step for many operations. Thus making it available to the user permits more rapid computation when several operations such as solving equations and determinant evaluation are to be performed using the same matrix. However, modern hardware is so fast that this would only seem to be necessary for very large sets of equations and these are not the target of these simple facilities. Moreover, adding explicit LU decomposition introduces complexity for the user.

Consideration was also given to a fuller implementation of the BLAS. However, this seems out of place in a language standard since it would be extremely long and specialized. Such a fuller interface to the BLAS could be provided in additional child packages. The goal here has been to provide convenient access to simple type declarations and the more commonly required operations that are not trivial for the user to program. These operations can of course be implemented as a binding to an implementation of part of the BLAS).

generic

   type Real is digits <>;

package Ada.Numerics.Generic_Real_Arrays is

pragma Pure (Generic_Real_Arrays);

type Real_Vector is array (Integer range <>) of Real'Base;

type Real_Matrix is array (Integer range <>, Integer range <>) of Real'Base;

function "+"   (Right : Real_Vector) return Real_Vector;

function "-"   (Right : Real_Vector) return Real_Vector;

function "abs" (Right : Real_Vector) return Real_Vector;

function "+" (Left, Right : Real_Vector) return Real_Vector;

function "-" (Left, Right : Real_Vector) return Real_Vector;

function "*" (Left, Right : Real_Vector) return Real'Base;

function "*" (Left : Real'Base; Right : Real_Vector) return Real_Vector;

function "*" (Left : Real_Vector; Right : Real'Base) return Real_Vector;

function "/" (Left : Real_Vector; Right : Real'Base) return Real_Vector;

function Unit_Vector (Index : Integer;
                      Order : Positive;
                      First : Integer := 1) return Real_Vector;

function "+"   (Right : Real_Matrix) return Real_Matrix;

function "-"   (Right : Real_Matrix) return Real_Matrix;

function "abs" (Right : Real_Matrix) return Real_Matrix;

function Transpose (X : Real_Matrix) return Real_Matrix;

function "+" (Left, Right : Real_Matrix) return Real_Matrix;

function "-" (Left, Right : Real_Matrix) return Real_Matrix;

function "*" (Left, Right : Real_Matrix) return Real_Matrix;

function "*" (Left, Right : Real_Vector) return Real_Matrix;

function "*" (Left : Real_Vector; Right : Real_Matrix) return Real_Vector;

function "*" (Left : Real_Matrix; Right : Real_Vector) return Real_Vector;

function "*" (Left : Real'Base; Right : Real_Matrix) return Real_Matrix;

function "*" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix;

function "/" (Left : Real_Matrix; Right : Real'Base) return Real_Matrix;

function Solve (A : Real_Matrix; X : Real_Vector) return Real_Vector;

function Solve (A, X : Real_Matrix) return Real_Matrix;

function Inverse (A : Real_Matrix) return Real_Matrix;

function Determinant (A : Real_Matrix) return Real'Base;

function Identity_Matrix (Order            : Positive;
                          First_1, First_2 : Integer := 1)
   return Real_Matrix;

function Unit_Matrix (Order            : Positive;
                      First_1, First_2 : Integer := 1) return Real_Matrix;

function Eigenvalues (A : Real_Matrix) return Real_Vector;

procedure Eigensystem (A       :     Real_Matrix;
                       Values  : out Real_Vector;
                       Vectors : out Real_Matrix);

private

   --  Implementation-defined ...

end Ada.Numerics.Generic_Real_Arrays;