This is to announce a new version of a Linear Algebra and
Optimization package, LinAlg v4.3, and to exhibit various matrix
streams it supports. The new version adds SVD and FFT pieces that were
missing from LinAlg 4.1. A web page (below) gives an overview of the
package and its features; the README file talks about them in far
greater detail. This article will try to demonstrate only one paradigm
of the package: Matrix streams.
Matrix streams have been introduced in the previous version of
LinAlg, but it's version 4.3 that enhances them and expands their
scope and applications. I have kept in mind two goals, efficiency
_and_ safety, and tried to satisfy them both.
In the pecking order, Matrix streams are considered
views. There are two kinds of sequential views: ElementWise streams
are transient views that exist only for the duration of an elementwise
action. Hence these views don't require a "locking" of a matrix. On
the other hand, LAStream and the ilk are more permanent views. An
application traverses the stream at its own convenience, bookmarking
certain spots and possibly returning to them later. An LAStream does
assert a shared "lock" on the matrix, to prevent the grid of matrix
values from moving or disappearing.
LinAlg 4.1 permitted creation of a Matrix stream -- LAStream --
over an entire Matrix, or its separate row, column, or the
diagonal. It is possible now to make an LAStream over an arbitrary
rectangular block of a matrix (including the whole matrix and a single
matrix element, and all other block sizes in between). This makes it
trivial to assign a rectangular block of one matrix to a rectangular
(sub)block of another.
LAStreams are now 'seek'able, from their beginning, end, or the
current position. An offset can be a positive, zero, or a negative
number. On no occasion one may position a stream prior to its lower
boundary. Seeking past the upper boundary does not crash the program
but merely sets an EOF condition, which can later be checked. All
flavors of LAStreams are seeakable -- a unit stride stream, an
arbitrary stride stream, and a block stream. Here's one application of
a 'fast forwarding' of a stream, which shows off a 'triangular
stream':
// Filling in the upper triangle of the matrix
// m[i,j] = seed + (i-1) + j*(j-1)/2; j>= i
{
LAStreamOut m_str(m);
for(register int j=0, curr=seed; j<m.q_ncols(); j++)
{
for(register int i=0; i<=j; i++)
m_str.get() = curr++;
m_str.ignore(m.q_ncols()-j-1);
}
}
As with regular file streams, a Matrix stream can be either writable,
or read-only. One may construct a writable LAStream only over a
non-const Matrix or Matrix reference.
It's now possible to create a stream spanning over a (part) of
another stream: subrange a stream. A substream is always created
inlined, it does not allocate any heap storage, and is _safe_. That
is, a substream may never extend its parent. A substream does not care
if the parent stream dies, as a substream asserts its own shared lock
on the original container. If a parent stream is writable, a child may
be created as read-only (but not the other way around!). Substream
creation is rather flexible and expressive, as the following snippet
shows:
static void upper_to_lower_triangle_reflector2(Matrix& m)
{
// must be a square matrix
assert( m.q_ncols() == m.q_nrows() );
for(register int i=0; i<m.q_ncols()-1; i++)
{
LAStreamOut out_str(MatrixColumn(m,i+m.q_col_lwb()),IRange::from(i+1));
LAStrideStreamIn in_str(ConstMatrixRow(m,i+m.q_col_lwb()),
IRange::from(i+1));
while( !out_str.eof() )
out_str.get() = in_str.get();
}
}
Previously I have shown that matrix inverse and determinant evaluation
can be done only sequentially, that is, without random access to
matrix elements. Although matrix inverse can be of value by itself
(esp. in Linear Regression, to find a correlation matrix), it doesn't
make the great example of using streams. The new version shows that
Householder bidiagonalizations can also be performed with only serial
access to matrices A, U, V in question. Here's a snippet from a
SVD::left_householder() method (contrasted with a similar snippet from
the previous version of LinALg, which did use random access).
Before (_un_serialized algorithm)
for(j=i+1; j<=N; j++) // Transform i+1:N columns of A
{
MatrixColumn Aj(A,j);
REAL factor = 0;
for(k=i; k<=M; k++)
factor += UPi(k) * Aj(k); // Compute UPi' * A[,j]
factor /= beta;
for(k=i; k<=M; k++)
Aj(k) -= UPi(k) * factor;
}
Here's a serialized algorithm from LinAlg 4.3, which employs
substreams that span only a part of matrix' columns:
IRange range = IRange::from(i - A.q_col_lwb());
LAStreamOut UPi_str(MatrixColumn(A,i),range);
...
for(register int j=i+1; j<=N; j++) // Transform i+1:N columns of A
{
LAStreamOut Aj_str(MatrixColumn(A,j),range);
REAL factor = 0;
while( !UPi_str.eof() )
factor += UPi_str.get() * Aj_str.get(); // Compute UPi' * A[,j]
factor /= beta;
for(UPi_str.rewind(), Aj_str.rewind(); !UPi_str.eof(); )
Aj_str.get() -= UPi_str.get() * factor;
UPi_str.rewind();
}
Unlike LinAlg4.1, this version of the library includes both
SVD and FFT, which were somewhat modified to take advantage of the new
features. Although FFT was hardly touched at all, SVD code was
_slightly_ changed to completely _avoid_ any random access to matrices
A, U, V. This shows that indeed many linear algebra algorithms -- even
as advanced as SVD -- could be implemented with only a sequential,
stream-like access to matrices or their parts.
Downloads:
http://pobox.com/~oleg/ftp/packages/LinAlg.tar.gz
The new README file, which includes detailed discussions of the new
concepts, examples, and a change log, can also be viewed separately:
http://pobox.com/~oleg/ftp/LinAlg.README.txt
Also viewable separately are the two interface files, LinAlg.h and
LAStreams.h
http://pobox.com/~oleg/ftp/packages/LinAlg.h
http://pobox.com/~oleg/ftp/packages/LAStreams.h
The FFT code for historical reasons is separated into:
http://pobox.com/~oleg/ftp/packages/fft.tar.gz
It is dependent upon the LinAlg43 code.
Brief description and overview:
http://pobox.com/~oleg/ftp/packages.html
Platforms supported:
gcc 2.7.2.1 Linux 2.0.33 i586
egcs-1.1.1 (optimization off, no '-freg-struct-return' gcc flag) Linux i686;
gcc 2.8.1 (optimization off) on HP-PA, Solaris/Sparc;
SGI's native compiler (full optimization off);
Visual C++ 5.0 WinNT 4.0;
Metrowerk's CodeWarrior C++, BeOS Release 4.
While the package correctly compiles with g++.2.8.1 and
egcs-1.1.1, I strongly recommend against using this compiler:
unfortunately, the only way it can handle my code is with optimization
off. Otherwise, the compiler quickly crashes with internal compiler
errors. It was my experience in general that g++ 2.8.1 is a somewhat
flaky version (as was g++ 2.4.0): unable to perform uncomplicated type
inferences and often tripping in its own extensions. Some compilers
consider themselves too smart and emit a _number_ of warnings
concerning 'initialization of non-const ... while passing ...'
etc. when compiling LinAlg4.3. Please disregard them.
This archive was generated by hypermail 2b29 : Wed Feb 20 2002 - 03:20:07 EST