readme

basic linear algebra classes and applications (SVD,interpolation, multivariate optimization)

them easier to build. Many linear algebra operations actually require
only sequential access to matrix elements. Thus the column index and other
direct access facilities are often redundant.

	If one does need to access arbitrary elements of a matrix, the package
offers several choices. One may use a special MatrixDA view, which is
designed to provide an efficient direct access to a matrix grid. The MatrixDA
view allocates and builds a column index, and uses it to efficiently
compute a reference to the (i,j)-th element:

  Matrix ethc = m;
  MatrixDA eth(ethc);
  for(register int i=eth.q_row_lwb(); i<=eth.q_row_upb(); i++)
    for(register int j=eth.q_col_lwb(); j<=eth.q_col_upb(); j++)
       eth(i,j) *= v(j);

see the validation code vmatrix*.cc, vvector.cc for more details. If creating
of an index seems like overkill, one may use MatrixRow, matrixColumn or
MatrixDiag. For example, the following statements are all equivalent:
  MatrixRow(m,2)(3) = M_PI;
  MatrixColumn(m,3)(2) = M_PI;
  MatrixDA ma(m); ma(2,3) = M_PI;
  (MatrixDA(m))(2,3) = M_PI;

Unlike Matrix itself, a Vector and all matrix views (MatrixRow,
MatrixColumn, MatrixDiag) allow direct access to their elements.

	In most of the cases, MatrixDA behaves like a matrix. On occasions
when one needs to get hold of a real matrix, he can always
use MatrixDA::ref() or MatrixDA::get_container() methods:
  Matrix mc(-1,10,1,20);
  MatrixDA m(mc);
  m.get_container().clear();
  verify_element_value(m,0);

	Again, it is strongly encouraged to use stream-lined matrix
operations as much as possible. As an example, finding of a matrix inverse
and computing of a determinant are implemented in this package using only
sequential matrix access. It is possible indeed to avoid direct access.


7. Const and Writable views

	There are two flavors of every matrix view: a const and a writable
views, for example: ConstMatrixDA and MatrixDA, ConstMatrixColumn and
MatrixColumn, ElementWiseConst and ElementWise, LAStreamIn and LAStreamOut.
A const view can be constructed around a const matrix reference. A const
view provides methods that do not alter the matrix: methods that
query matrix elements, compare them, compute norms and other cumulative
quantities (sum of squares, max absolute value, etc). A writable view is
a subclass of the const view. That means that a writable view allows all
the query/comparison operations that the corresponding const view implements.
In addition, a writable view permits modification of matrix elements, via
assignment or other related operations (+=, *=, etc). Needless to say one
needs a non-const matrix reference to construct a writable view:

  cout << "   compare (m = pattern^2)/pattern with pattern\n";
  m = pattern; to_every(m1) = pattern;
  to_every(m).sqr();
  assert( of_every(m).max_abs() == pattern*pattern );
  to_every(m) /= of_every(m1);
  verify_element_value(m,pattern);

to_every() makes a writable ElementWise view, while of_every() makes a
ElementWiseConst view.


8. Never return non-trivial objects (matrices or vectors)

Danger: For example, when the following snippet 
  Matrix foo(const int n)
  { Matrix foom(n,n); fill_in(foom); return foom; }
  Matrix m = foo(5);
runs, it constructs matrix foo:foom, copies it onto stack as a return
value and destroys foo:foom. The return value (a matrix) from foo() is
then copied over to m via a copy constructor. After that, the return value
is destroyed. The matrix constructor is called 3 times and the
destructor 2 times. For big matrices, the cost of multiple
constructing/copying/destroying of objects may be very large. *Some*
optimized compilers can do a little better. That still leaves at least
two calls to the Matrix constructor. LazyMatrices (see below) can construct
a Matrix m "inplace", with only a _single_ call to the constructor.


9. Lazy matrices

	Instead of returning an object return a "recipe" how
to make it. The full matrix would be rolled out only when and where
it is needed:
  Matrix haar = haar_matrix(5);
haar_matrix() is a *class*, not a simple function. However similar
this looks to returning of an object (see the note above), it is
dramatically different. haar_matrix() constructs a LazyMatrix, an
object just of a few bytes long. A special "Matrix(const LazyMatrix&
recipe)" constructor follows the recipe and makes the matrix haar()
right in place. No matrix element is moved whatsoever!

Another example of matrix promises is
  REAL array[] = {0,-4,0,  1,0,0,  0,0,9 };
  test_svd_expansion(MakeMatrix(3,3,array,sizeof(array)/sizeof(array[0])));

Here, MakeMatrix is a LazyMatrix that promises to deliver a matrix
filled in from a specified array. Function test_svd_expansion(const Matrix&)
forces the promise: the compiler makes a temporary matrix, fills
it in using LazyMatrix's recipe, and passes it out to test_svd_expansion().
Once the function returns, the temporary matrix is disposed of.
All this goes behind the scenes. See vsvd.cc for more details (this is
where the fragment was snipped from).

	Lazy matrices turned out so general and convenient that they subsumed
glorified constructors (existed in the previous version of the LinAlg). It is
possible to write now:
  Matrix A = zero(B);
  Matrix C = unit(B);
  Matrix D = transposed(B);
  A = unit(B);
  D = inverse(B);
In all these cases, the result is computed right in-place, no temporary
matrices are created or copied.

	Here's a more advanced example, 
  Matrix haar = haar_matrix(5);
  Matrix unitm = unit(haar);
  Matrix haar_t = transposed(haar);
  Matrix hht = haar * haar_t;		// NB! And it's efficient!
  Matrix hht1 = haar; hht1 *= haar_t;
  Matrix hht2 = zero(haar); hht2.mult(haar,haar_t);
  verify_matrix_identity(unitm,hht);
  verify_matrix_identity(unitm,hht1);
  verify_matrix_identity(unitm,hht2);

	With lazy matrices, it is possible to express a matrix multiplication
in a "natural way", without the usual penalties:
  C = A * B;
Here again, the star operator does not actually multiply the matrices. It
merely constructs a lazy matrix, a promise to compute the product (when the
destination of the product becomes known). The promise is then forced by the
assignment to matrix C. A check is made that the dimensions of the promise
are compatible with those of the target. The result will be computed right
into matrix's C grid, no temporary storage is ever allocated/used.

	As even more interesting example, consider the following snippet
from the verification code:
  verify_matrix_identity(m * morig,unit(m));
Here both argument expressions compute matrix promises; the promises are
forced by parameter passing, so that the function can compare the resulting
matrices and verify that they are identical. The forced matrices are 
discarded afterwards (as we don't need them). Note how this statement leaves
no garbage.

	A file sample_adv.cc tries to prove that LazyMatrices do indeed
improve performance. You can see this for yourself by compiling and
running this code on your platform. LazyMatrices are typically 2 to
three (Metrowerks C++, BeOS) times as efficient.


10. Lazy matrices and expression templates

        Lazy matrices and expression templates are a bit different.
Expression templates are a _parser_ of expressions, an embedded language.
They parse an expression like
        C = A*B;
at compile time and generate the most efficient, inlined code. This
however places an immense burden on a compiler. For one thing, the
compiler has to optimize out _very_ many temporary objects of various types,
which emulate Expression templates parser's stack. The compiler should
be able to automatically instantiate templates; gcc until recently couldn't,
and even now I won't dare to ask it to. Furthermore, one has to be very
careful writing expression templates: every type of expression
(multiplication, addition, parenthesized expression, left-right
multiplication by a scalar (a double, float,...)), etc. requires 
its own template. One has to make sure that multiplication
is indeed performed ahead of addition, according to usual operation
priorities, as in D = A + B*C; Templates are not the best tool to write
parsers to start with. The very need to write one's own parser in C++ (which
already has a parser) appears rather contrived: this may show that C++
might not be the right language for such problems after all.

        Lazy matrices do not require any templates at all; multiplication
itself as in
        C = A*B;
is performed by a separate function, which may be a library function
and may not even be written in C++: as a typical BLAS, matrix multiplication
can be coded in assembly to take the best advantage of certain architectures.

        Lazy matrices is a technique of dealing with an operation like
A*B when its target is not yet known. In the standard C++ paradigm,
the expression is evaluated and the result is placed into a temporary.
Lazy matrices offer a way to delay the execution until the target,
the true and the final recipient of the product, becomes known.
One doesn't need any matrix temporary then. Lazy computation is the
standard fare of functional programming; in Haskell, all computation is
lazy. That is, a Haskell code isn't in a hurry to execute an operation,
it waits to see if the results are really necessary, and if so,
where to place them. It's interesting  that a similar technique is
possible in C++ as well.

        As lazy matrices and expression templates are rather independent,
one can mix and match them freely.


11. Accessing row/col/diagonal of a matrix without much fuss
(and without moving a lot of stuff around)

  Matrix m(n,n); Vector v(n); 
  to_every(MatrixDiag(m)) += 4;	// increments the diag elements of m
				// by 4

  v = of_every(ConstMatrixRow(m,1));	// Saves a matrix row into a vector

  MatrixColumn(m,1)(2) = 3; // the same as m(2,1)=3, but does not
			    // require constructing of MatrixDA

  (MatrixDiag(m))(3) = M_PI; // assigning pi to the 3d diag element

Note, constructing of, say, MatrixDiag does *not* involve any copying
of any elements of the source matrix. No heap storage is used either.


12. Nested functions

	For example, creating of a Hilbert matrix can be done as follows:

  void foo(const Matrix& m)
  {
    Matrix m1 = zero(m);
    struct MakeHilbert : public ElementAction
    {
      void operation(REAL& element, const int i, const int j)
	{ element = 1./(i+j-1); }
    };
    m1.apply(MakeHilbert());
  }

A special method Matrix::hilbert_matrix() is still more optimal, but not
by the whole lot. The class MakeHilbert is declared *within* a function
and is local to that function. This means one can define another MakeHilbert
class (within another function or outside of any function - in the global
scope), and it will still be OK.

Another example is application of a simple function to each matrix element

 void foo(Matrix& m, Matrix& m1)
 {
  class ApplyFunction : public ElementWiseAction
  {
    double (*func)(const double x);
    void  operator () (REAL& element) { element = func(element); }
    public: ApplyFunction(double (*_func)(const double x)) : func(_func) {}
  };
  to_every(m).apply(ApplyFunction(sin));
  to_every(m1).apply(ApplyFunction(cos));
 }

Validation code vmatrix.cc and vvector.cc contains a few more examples
of this kind (especially vmatrix1.cc:test_special_creation())


13. SVD decomposition and its applications

Class SVD performs a  Singular Value Decomposition of a rectangular matrix
A = U * Sig * V'. Here, matrices U and V are orthogonal; matrix Sig is a
diagonal matrix: its diagonal elements, which are all non-negative, are
singular values (numbers) of the original matrix A. In another interpretation,
the singular values are eigenvalues of matrix A'A.

Application of SVD: _regularized_ solution of a set of simultaneous
linear equations Ax=B.  Matrix A does _not_ have to be a square matrix.
If A is a MxN rectangular matrix with M>N, the set Ax=b is obviously
overspecified. The solution x produced by SVD_inv_mult would then be  
the least-norm solution, that is, a least-squares solution.
Note, B can be either a vector (Mx1-matrix), or a "thick" matrix. If B is
chosen to be a unit matrix, then x is A's inverse, or a pseudo-inverse if
A is non-square and/or singular.
An example of using SVD:
	SVD svd(A);
	cout << "condition number of matrix A " << svd.q_cond_number();
	Vector x = SVD_inv_mult(svd,b);		// Solution of Ax=b
Note that SVD_inv_mult is a LazyMatrix. That is, the actual computation
occurs not when the object SVD_inv_mult is constructed, but when it's
required (in an assignment).


14. Functor as a function pointer

Package's fminbr(), zeroin(), hjmin() are functions of "the higher
order": they take a function, any function, and try to find out a
particular property of it - a root or a minimum. In LinAlg 3.2 and
earlier, this function under investigation was specified as a pointer
to a C/C++ function. For example,

static int counter;
static double f1(const double x)		// Test from the Forsythe book
{
  counter++;
  return (sqr(x)-2)*x - 5;
}
main()
{
  counter = 0;
  const double a = 0, const double b = 1.0;
  const double found_location_of_min = fminbr(a,b,f1);
  printf("\nIt took %d iterations\n",counter);
}

Note that function f1(x) under investigation had a side-effect: it
alters a 'counter', to count the total number of function's
invocations. Since f1(x) is called (by fminbr()) with only single
argument, the 'counter' has to be declared an external variable. This
makes it exposed to other functions of the package, which may
inadvertently change its value. Also, the main() function should know
that f1() uses this parameter 'counter', which main() must initialize
properly (because the function f1() obviously can't do that itself)

	The new, 4.0 version of the package provides a more general
way to specify a function to operate upon. A math_num.h file declares
a generic function class, a functor, which takes a single
floating-point argument and returns a single FP (double) number as the
result:

struct UnivariateFunctor {
  virtual double operator() (const double x) = 0;
};

The user should specialize this abstract base class, inheriting from
it and specifying "double operator()" to be a function being
investigated.  The user may also add his own data members and methods
to his derived class, which would serve as an "environment" for the
function being optimized. The user should then instantiate this class
(thus initializing the environment), and pass the instance to
fminbr().  The environment will persist even after fminbr() returns,
so the user may use accumulated results it may contain.  For example,

		// Declaration of a particular functor class
class F1 : public UnivariateFunctor
{
protected:
  int counter;                  // Counts calls to the function
  const char * const title;
public:
  ATestFunction(const char _title []) : counter(0), title(_title) {}
  double operator() (const double x) { counter++; return (sqr(x)-2)*x - 5; }
  int get_counter(void) const	{ return counter; }
};

main()