GSoC 2012 Proposal

cristicbz edited this page Mar 28, 2012 · 4 revisions
Clone this wiki locally

Continued Work on a D Linear Algebra library (SciD - std.linalg)

This is a proposal to continue my work on a D linear algebra library, a task which I started last year, mentored by David Simcha and co-mentored by Fawzi Mohamed and Andrei Alexandrescu.

The main aim of this proposal would be to bring my fork of SciD to a version in which it could be considered for inclusion in Phobos as std.linalg (or some possibly better name). This means that unlike last year, the effort would be much more focus on getting a smaller set of features working with a high degree of reliability and possibly moving the more advanced features into an 'experimental' add-on to the library, if not remove them altogether.

To this effect, I propose a redesign of the library's interface inspired by the C++ library Eigen keeping the core code unchanged and adding support for fixed-size matrices, value arrays (matrices that support less features, but allow arbitrary element types), creating a generic model for decompositions (which, of course, wrap LAPACK functions). As mentioned before, I also believe that the removal (or at least deprecation) of certain features would be beneficial - views in particular.

A suggestion made by David Simcha, with which I entirely agree, is to set a soft deadline at the end of July after which no more features would be added, spending August getting the library ready for inclusion in Phobos.

Current Situation and Justification for changes

The aims of my last GSoC project last year expanded significantly during last summer, from a contribution to the existing SciD library to an essentially new library, mostly written from scratch, supporting many more features than originally planned. A significant proportion of the time was also spent isolating, working around and, occasionally, fixing compiler bugs, and soon my mentor and I realised that our goals had become too ambitious for the allotted time. Many goals were achieved - expressive matrix & vector types, expression template construction and evaluation using BLAS/LAPACK, naive implementations for unavailable libraries etc. - but many were not and at the end of the summer the project was left in a rough, alpha state. I had hoped to be able to complete the work during the year, but unfortunately I simply did not have as much spare time as I would have wished. The foundation created last summer, however, is ideal for building upon; the core of the library: memory management, expression evaluation, naive implementations are thoroughly tested and perfectly usable.

Furthermore, during the past year I've done much more numerical programming than I have ever done before, using MATLAB and Eigen extensively. The experience I gained puts me in a much better position to estimate the priorities of different features and I hope it will allow me to better judge a good library design.

Proposed Changes and Additions

This can also be read as a roadmap since I plan on making these changes more-or-less in order:

  1. Remove views and packed storage temporarily. In short, this is because they add a lot of complexity to the system for different reasons. I believe a better design would be to build them on top of other concepts rather than as an integral part of the library. See the relevant section (Views and Packed Storage) below for justification.
  2. Change the matrix & vector types, adding fixed-sized matrix support in the process. The original interface had some good ideas and some not as good. I would particularly like to remove the difference between Vector and Matrix classes, trying to avoid duplication between their implementation was pointlessly difficult. The proposed design also adds fixed-sized matrices as a fundamental part of the library. Example:
Matrix!( double, 2, 2 )                    twoByTwoMatrix;
Matrix!( float, DynamicSize, DynamicSize ) dynamicMatrix;
Matrix!( real, 1, DynamicSize )            realRowVector;
Matrix!( Complex!double, DynamicSize, 1 )  complexColumnVector;

alias Matrix!( double, 2,           1           ) Vector2d;
alias Matrix!( double, 3,           1           ) Vector3d;
alias Matrix!( double, 4,           1           ) Vector4d;
alias Matrix!( double, DynamicSize, 1           ) VectorXd;
alias Matrix!( double, DynamicSize, DynamicSize ) MatrixXd;

Matrices are by default column-major, for row-major matrices the fourth template parameter can be set to StorageOrder.RowMajor. Example: Matrix!( double, DynamicSize, DynamicSize, StorageOrder.RowMajor ). For convenience, there should be a type RowMajorMatrix!double which is equivalent to the latter. 3. Add value arrays (or numeric arrays, we can come up with a good name). These are light-weight matrices which support element-wise operations, arbitrary element types, higher dimensionality (greater than two). Conversion to and from matrices is free, performed with the .array() and .matrix() members correspondingly. The array/matrix returned shares memory with the original in a copy-on-write fashion. Example:

Array!( double, 2, 2 ) twoByTwoArray;
twoByTwoMatrix.set( 1, 2,
                    3, 4 );
writeln( "Element-wise multiplication: ", twoByTwoArray * twoByTwoArray );
writeln( "Matrix multiplication: ", twoByTwoArray.matrix() * twoByTwoArray.matrix() );

twoByTwoMatrix += 1;
writeln( "Add 1 to all elements: ", twoByTwoMatrix );
  1. Add reductions, partial reductions, and broadcasting for matrices and arrays. Reductions refer to operations such as computing the sum/product/min/max of elements. Partial reductions are column/row-wise reduction (computing a row vector of the sum along each column of the matrix). Broadcasting is the complement of partial reductions - it corresponds to MATLAB's repmat() - the best example is adding a row vector to all the rows of a matrix. Examples:
auto mat = MatrixXd( 3, 3 );
mat.set( 1, 2, 3,
         4, 5, 6,
         7, 8, 9 );

// Reductions
writeln( "Sum and product: ", mat.sum(), ", ", );

auto indexValue = mat.argmax();
writeln( "Argmax: mat[", indexValue.index( 0 ), ", ", indexValue.index( 1 ), "] = ",  indexValue.value ); 

// Partial Reduction (outputs: "Sum along columns and along rows: [12 15 18], [6; 9; 12]")
writeln( "Sum along columns and along rows: ", mat.rowwise().sum(), ", ", mat.colwise().sum() );

// Brodcasting (.t performs transposition which is neccessary, since vec is a column vector)
// Outputs: "Add [1; 2; 3] to all rows: [2 4 6 ; 5 7 9 ; 8 10 12]"
auto vec = VectorXd( 3 );
vec.set( 1,2,3 );
writeln( "Add ", vec, " to all rows: ", mat.rowwise() + vec.t );

The .rowwise() and .colwise() functions can be defined elegantly as returning a value array (see 3) of rows and columns respectively - this automatically results in the desired functionality.

  • Note: Calling .rowwise() or .colwise() is essentially free. The return type of Matrix!( double, nRows, nColumns ).rowwise() is Array!( Matrix!( double, 1, nColumns ), nRows, 1 ). The returned array shares memory with the original memory in the usual copy-on-write way and therefore there is no memory allocation: pointer assignment, int-increment and int-decrement is all that happens.
  1. Add support for boolean arrays. These are the result of comparison and equality operators on matrices or arrays, as well as logical operators on other boolean arrays. They can be reduced using .any(), .all(), .sum() (for counting) e.g. (mat > 3).all() returns true if all elements in mat are greater than three.
  • Note: Again, as in the case of broadcasting, boolean arrays are free. This is because no operations are evaluated until assignment - including logical operations. So, writing bmat = (mat > 3) would allocate the memory necessary to save the matrix, while (mat > 3).all() would be translated into a loop equivalent to:
foreach( auto column ; mat )
    foreach( auto elem ; column )
        if( !(elem > 3) )
            return false;
return true;
  1. Add support for interoperation with D built-in arrays (or pointers). This will be done through the External!( MappedType, StrideType=DefaultStride ) type. For example, to map a double[9] to a Matrix!( double, 3, 3 ):
double[ 9 ] arr = [ 1, 2, 3, 4, 5, 6, 7, 8, 9 ]
auto ext1 = External!( Matrix!( double, 3, 3 ) )( arr );
auto ext2 = externalMatrix!( 3, 3 )( arr ); // equiv. to ext1, convenience function
auto ext3 = externalMatrix( arr, 3, 3 );    // using dynamic size rather than compile-time
ext1[ 0, 0 ] = 42;
assert( arr[ 0 ] == 42 && ext1[ 0, 0 ] == 42 && ext2[ 0, 0 ] == 42 && ext3[ 0, 0 ] == 42 );

StrideType specifies the distance in memory between matrix elements and should be a specialisation of the template Stride( size_t InnerStride = DynamicSize, size_t OuterStride = DynamicSize ). The inner stride is the pointer increment between two consecutive elements in a column for column-major matrices or in a row for a row-major matrices. The outer stride is the pointer increment between two columns in column-major matrices or between two rows for row-major matrices. The Stride type, similar to matrices, allows for compile-time or dynamic strides - if the template parameters are set to DynamicSize, then the stride values can be set by passing a Stride object to the External constructor:

double[ 9 ] arr = [ 1, 2, 3,
                    4, 5, 6,
                    7, 8, 9 ];
auto v1 = External!( Matrix!( double, 1, 3 ), Stride!(3, 1) )( arr ); // first column of arr [1, 4, 7]
auto v2 = externalMatrix!( 1, 3, Stride!( 3, 1 ) )( arr );            // equivalent to v2
auto v3 = externalMatrix!( 1, 3 )( arr, Stride(3, 1) );               // using dynamic stride
auto v4 = externalMatrix( arr, 1, 3, innerStride(3) );                // using dynamic stride and size
                                                                      //   innerStride( 3 ) == Stride(3, 1)
  1. Add a general model for decompositions. Through a decomposition one can compute multiple properties, so I propose a version where decompositions are represented by types conforming to the Decomposition concept:
template isDecomposition( Decomposition ) {
    enum isDecomposition =
        __traits( compiles, (){
            auto d   = Decomposition( MatrixXd(2,2) );
            auto det = d.determinant();
            auto inv = d.inverse();
            auto sol = d.solve( VectorXd( 2 ) );

Then to use a decomposition one can use this kind of syntax:

Cholesky!MatrixXd ch = cholesky( mat ); // return type specified for clarity
writeln( ch.inverse() );
writeln( ch.solve( vec ) );

// Of course, one does not need to save the decomposition to use the result:
writeln( lu( mat ).inverse() );

// As in the current version of the library, the following are evaluated to the same LAPACK calls:
writeln( inv( mat ) * vec );
writeln( lu( mat ).solve( vec ) );

Removed Features

Packed Storage

Right now, matrices wrap four level of abstraction, which allows 'customization' every step of the way. Except it doesn't. In an effort to write DRY, yet customisable code, serious sacrifices were made in terms of usability. Customising any of the levels is a pain, requiring intimate knowledge of multiple concepts (Storage, Container, Operand and Promote, in particular). This customisability was added mostly to allow implementation of packed storage types, diagonal matrices and the like with minimal code repetition.

In retrospect, this sacrifice was not worth the results. The new design should implement packed storage types as different types altogether, keeping the dense, general matrices self-contained. Operations between general matrices and packed matrices should be handled by the latter (opBinaryRight etc.). In all, I believe I overestimated the importance of packed storage.


Why are views evil? Well, I didn't say they were evil, but since you asked: a view in SciD is a reference type that refers to a matrix/vector, matrix/vector slice or matrix diagonal. The interesting bit is that views share ownership of memory between each other as well as with matrix/vector objects - you cannot get a dangling view, even if the last Matrix/Vector object to which the view refers is destroyed, since the view owns the memory. Add to this the fact that Matrix/Vector types are copy-on-write, and you end up needing two levels of reference counting (see the wiki page at, for a more in-depth explanation). This again causes over-complication of the design and as a side effect, matrix assignments have to allocate memory on the heap.

About Me

I am a 3rd year student at the University of Edinburgh, UK where I am double-majoring in Artificial Intelligence and Computer Science where I have an average mark of 92% so far. While C++ is my native language and while I am quite proficient in Haskell, Python and Java, I have always had soft spot for D and tried to write any personal scrap of code in it. As mentioned, I worked on SciD last summer, building the library that exists now - this I hope provides some notion on my ability with meta-programming and memory management in D, as well as my knowledge of BLAS and LAPACK.

Many of the courses I am taking (computer vision, machine learning and natural language processing courses, in particular) are numerics-heavy and provided me with some insight into which features are important, what priorities to set and what trade-offs to make in a numerics package. This also allowed me to gain more experience working with existing solutions - MATLAB and Eigen to be precise.