This project has moved and is read-only. For the latest updates, please go here.

Symmetric Dense Matrices

Mar 9, 2011 at 4:46 PM
Edited Mar 9, 2011 at 4:57 PM

I am working on a group of classes that handle Symmetric dense matrices. I need some input on structure/architecture/design/whatever_this_is_called to make it fit nicely.

Symmetric matrices provide an immediate memory advantage. A (square) dense matrices has n^2 values whereas a symmetric dense needs n*(n+1)/2 (roughly half+) Also, we save a lot on calculations (eg when creating it we only need to set half+ the values etc). It would be very advantageous to be able to specify when a matrix is symmetric. By convention, I am saving the upper triangle of the matrix row-wise (the partial row from the diagonal element and to the right, then the next partial row etc).

A symmetric matrix can be saved with 2 major ways:
1) Vector form. Since we know exactly how many elements are saved on each row, it is easy to calculate where an element [i,j] is located on the vector and access it.
2) Jagged Array form(eg double[][]). Each sub-array contains exactly as many elements as needed. Eg the first one contains n elements, the second one contains n-1 elements etc.

I am strongly in favor of the second form because its a lot clearer and maintainable. However, I noticed that even for DenseMatrices MathNet uses double[]. I would think that having to perform (column * RowCount) + row to access every element would be less efficient than double[,] (even though [,] would make that sort of calculation internally) but I am sure this has been researched (I think Iridium used to use [,] or [][]). Please clarify which form is preferred.


Next issue: For the DenseMatrix class of Double/Single/Complex/Complex32 we got:
DenseMatrix : Matrix : Matrix<T>

Symmetric is a property of the matrix, so an interface ISymmetric seems natural. But since I cannot provide a body to the methods, I thought of creating a new abstract class next to Matrix. Currently I am thinking of:

SymmetricDenseMatrix : SymmetricMatrix : Matrix<T>

for each one of Double/Single/Complex/Complex32. This is would help in creating a: SymmetricSparseMatrix : SymmetricMatrix : Matrix<T>

How does this seem? I also need to take into account that the matrix is always square. Do I need a SquareMatrix<T>?

 

All input is greatly appreciated :).

Mar 9, 2011 at 6:06 PM
Edited Mar 9, 2011 at 6:07 PM

A symmetric matrix would be great addition.

>However, I noticed that even for DenseMatrices MathNet uses double[]. I would think that having to perform (column * RowCount) + row
>to access every element would be less efficient than double[,] 

We are trying to keep the matrix storage compatible with native BLAS/LAPACK matrix storage schemes. For dense matrices, the matrices are stored as a single dimensional array with LAPACK requiring the data to be stored in a column-major format (row-major might be an option now with the new C API for LAPACK, but I haven't checked and I don't know how well that is supported with different LAPACK implementations). I believe [,] is stored in contiguous memory in .NET so it can passed as a single dimensional array to the native libraries, but it will be in row-major form. Also in earlier versions of the .NET runtime, [,] based matrices were exceptionally slow. I don't know if the newer versions of the runtime have improved performance. It might be time to redo our benchmark comparing [],[][], and [,] ([] were the fastest and [,] were the slowest in the benchmarks we did years ago for dnAnalytics even with the [] indexing).  

As for the dense symmetric matrix storage, I think we need to use packed symmetric storage in a 1-D array to be compatible with the native libraries. I'll double check tomorrow, as well see if the CSB sparse storage scheme is compatible.

> Do I need a SquareMatrix<T>?
I wouldn't add one. 

Mar 10, 2011 at 3:17 PM

If we want to take advantage of native BLAS/LAPACK routines for the symmetric matrix, we should use packed storage. See this link for a description: http://software.intel.com/sites/products/documentation/hpc/mkl/webhelp/appendices/mkl_appB_MA.html

We are planning on using native routines for the sparse matrix, so we need to select a storage scheme that is compatible with MKL (and hopefully other providers). Here is a list of the MKL sparse storage formats: http://software.intel.com/sites/products/documentation/hpc/mkl/webhelp/appendices/mkl_appA_SMSF.html

We'd have to add a second sparse matrix type if we wanted to support CSB storage. User's could then select either based on their needs.

Mar 12, 2011 at 12:57 AM

Ok, I will change the symmetric storage to the scheme used in MKL.

Mar 24, 2011 at 8:06 PM

Are there methods in the LinearAlgebraProvider(s) for symmetric matrices?

Mar 25, 2011 at 6:00 AM

No. We are focusing on dense matrices for the first release. Fell free to modify the LinearAlgebraProvider interface to handle them.