newmatb.txt

科学和工程计算中使用统计功能开发工具包

//$$ newmatb.txt                                   Design


Copyright (C) 1991,2: R B Davies

Please treat this as an academic publication. You can use the ideas but
please acknowledge in any publication.



                      Notes on the design of the package
                      ==================================

Contents
========

What this is package for
What size of matrices?
Allow matrix expressions?
Which matrix types?
What element types?
Naming convention
Row and Column index ranges
Structure of matrix objects
Data storage - one block or several
Data storage - by row or by column or other
Storage of symmetric matrices
Element access - method and checking
Use iterators?
Memory management - reference counting or status variable?
Evaluation of expressions - use two stage method?
How to overcome an explosion in number of operations
Destruction of temporaries
A calculus of matrix types
Error handling
Sparse matrices
Complex matrices


I describe some of the ideas behind this package, some of the decisions
that I needed to make and give some details about the way it works. You
don't need to read this file in order to use the package.

I don't think it is obvious what is the best way of going about
structuring a matrix package. And I don't think you can figure this
out with "thought" experiments. Different people have to try out
different approaches. And someone else may have to figure out which is
best. Or, more likely, the ultimate packages will lift some ideas from
each of a variety of trial packages. So, I don't claim my package is an
"ultimate" package, but simply a trial of a number of ideas.

But I do hope it will meet the immediate requirements of some people who
need to carry out matrix calculations.


What this is package for
------------------------

The package is used for the manipulation of matrices, including the
standard operations such as multiplication as understood by numerical
analysts, engineers and mathematicians. A matrix is a two dimensional
array of numbers. However, very special operations such as matrix
multiplication are defined specifically for matrices. This means that a
"matrix" package tends to be different from a general "array" package.

I see a matrix package as providing the following

1.   Evaluation of matrix expressions in a form familiar to
     scientists and engineers. For example  A = B * (C + D.t());
2.   Access to the elements of a matrix;
3.   Access to submatrices;
4.   Common elementary matrix functions such as determinant and trace;
5.   A platform for developing advanced matrix functions such as eigen-
     value analysis;
6.   Good efficiency and storage management;
7.   Graceful exit from errors (I don't provide this yet).

It may also provide

8.   A variety of types of elements (eg real and complex);
9.   A variety of types of matrices (eg rectangular and symmetric).

In contrast an array package should provide

1'.  Arrays can be copied with a single instruction; may have some
     arithmetic operations, say + - and scalar + - * /, it may provide
     matrix multiplication as a function;
2'.  High speed access to elements directly and with iterators;
3'.  Access to subarrays and rows (and columns?);
6'.  Good efficiency and storage management;
7'.  Graceful exit from errors;
8'.  A variety of types of elements (eg real and complex);
9'.  One, two and three dimensional arrays, at least, with starting
     points of the indices set by user; may have arrays with ragged
     rows.

It may be possible to amalgamate these two sets of requirements to some
extent. However my package is definitely oriented towards the first set.

Even within the bounds set by the first set of requirements there is a
substantial opportunity for variation between what different matrix
packages might provide.

It is not possible to build a matrix package that will meet everyone's
requirements. In many cases if you put in one facility, you impose
overheads on everyone using the package. This both is storage required
for the program and in efficiency. Likewise a package that is optimised
towards handling large matrices is likely to become less efficient for
very small matrices where the administration time for the matrix may
become significant compared with the time to carry out the operations.

It is better to provide a variety of packages (hopefully compatible) so
that most users can find one that meets their requirements. This package
is intended to be one of these packages; but not all of them.

Since my background is in statistical methods, this package is oriented
towards the kinds things you need for statistical analyses.


What size of matrices?
----------------------

A matrix package may target small matrices (say 3 x 3), or medium sized
matrices, or very large matrices. A package targeting very small
matrices will seek to minimise administration. A package for medium
sized or very large matrices can spend more time on administration in
order to conserve space or optimise the evaluation of expressions. A
package for very large matrices will need to pay special attention to
storage and numerical properties.

This package is designed for medium sized matrices. This means it is
worth introducing some optimisations, but I don't have to worry about
the 64K limit that some compilers impose.


Allow matrix expressions?
-------------------------

I want to be able to write matrix expressions the way I would on paper.
So if I want to multiply two matrices and then add the transpose of a
third one I can write something like

   X = A * B + C.t();

I want this expression to be evaluated with close to the same efficiency
as a hand-coded version. This is not so much of a problem with
expressions including a multiply since the multiply will dominate the
time. However, it is not so easy to achieve with expressions with just +
and - .

A second requirement is that temporary matrices generated during the
evaluation of an expression are destroyed as quickly as possible.

A desirable feature is that a certain amount of "intelligence" be
displayed in the evaluation of an expression. For example, in the
expression

   X = A.i() * B;

where i() denotes inverse, it would be desirable if the inverse wasn't
explicitly calculated.


Which matrix types?
-------------------

As well as the usual rectangular matrices, matrices occuring repeatedly
in numerical calculations are upper and lower triangular matrices,
symmetric matrices and diagonal matrices. This is particularly the case
in calculations involving least squares and eigenvalue calculations. So
as a first stage these were the types I decided to include.

It is also necessary to have types row vector and column vector. In a
"matrix" package, in contrast to an "array" package, it is necessary to
have both these types since they behave differently in matrix
expressions. The vector types can be derived for the rectangular matrix
type, so having them does not greatly increase the complexity of the
package.

The problem with having several matrix types is the number of versions
of the binary operators one needs. If one has 5 distinct matrix types
then a simple package will need 25 versions of each of the binary
operators. In fact, we can evade this problem, but at the cost of some
complexity.


What element types?
-------------------

Ideally we would allow element types double, float, complex and int, at
least. It would be reasonably easy, using templates or equivalent, to
provide a package which could handle a variety of element types.
However, as soon as one starts implementing the binary operators between
matrices with different element types, again one gets an explosion in
the number of operations one needs to consider. Hence I decided to
implement only one element type. But the user can decide whether this is
float or double. The package assumes elements are of type Real. The user
typedefs Real to float or double.

In retrospect, it would not be too hard to include matrices with complex
matrix type.

It might also be worth including symmetric and triangular matrices with
extra precision elements (double or long double) to be used for storage
only and with a minimum of operations defined. These would be used for
accumulating the results of sums of squares and product matrices or
multistage Householder triangularisations.


Naming convention
-----------------

How are classes and public member functions to be named? As a general
rule I have spelt identifiers out in full with individual words being
capitalised. For example "UpperTriangularMatrix". If you don't like this
you can #define or typedef shorter names. This convention means you can
select an abbreviation scheme that makes sense to you.

The convention causes problems for Glockenspiel C++ on a PC feeding into
Microsoft C. The names Glockenspiel generates exceed the the 32
characters recognised by Microsoft C and ambiguities result. So it is
necessary to #define shorter names.

Exceptions to the general rule are the functions for transpose and
inverse. To make matrix expressions more like the corresponding
mathematical formulae, I have used the single letter abbreviations, t()
and i() .


Row and Column index ranges
---------------------------

In mathematical work matrix subscripts usually start at one. In C, array
subscripts start at zero. In Fortran, they start at one. Possibilities
for this package were to make them start at 0 or 1 or be arbitrary.
Alternatively one could specify an "index set" for indexing the rows and
columns of a matrix. One would be able to add or multiply matrices only
if the appropriate row and column index sets were identical.

In fact, I adopted the simpler convention of making the rows and columns
of a matrix be indexed by an integer starting at one, following the
traditional convention. In an earlier version of the package I had them
starting at zero, but even I was getting mixed up when trying to use
this earlier package. So I reverted to the more usual notation.


Structure of matrix objects
---------------------------

Each matrix object contains the basic information such as the number of
rows and columns and a status variable plus a pointer to the data
array which is on the heap.


Data storage - one block or several
-----------------------------------

In this package the elements of the matrix are stored as a single array.
Alternatives would have been to store each row as a separate array or a
set of adjacent rows as a separate array. The present solution
simplifies the program but limits the size of matrices in systems that
have a 64k byte (or other) limit on the size of arrays. The large arrays
may also cause problems for memory management in smaller machines.


Data storage - by row or by column or other
-------------------------------------------

In Fortran two dimensional arrays are stored by column. In most other
systems they are stored by row. I have followed this later convention.
This makes it easier to interface with other packages written in C but
harder to interface with those written in Fortran.

An alternative would be to store the elements by mid-sized rectangular
blocks. This might impose less strain on memory management when one
needs to access both rows and columns.


Storage of symmetric matrices
-----------------------------

Symmetric matrices are stored as lower triangular matrices. The decision
was pretty arbitrary, but it does slightly simplify the Cholesky
decomposition program.


Element access - method and checking
------------------------------------

We want to be able to use the notation A(i,j) to specify the (i,j)-th
element of a matrix. This is the way mathematicians expect to address
the elements of matrices. I consider the notation A[i][j] totally alien.
However I may include it in a later version to help people converting
from C. There are two ways of working out the address of A(i,j). One is
using a "dope" vector which contains the first address of each row.