This project has moved. For the latest updates, please go here.

Multiplication of different matrix-types

May 10, 2009 at 3:32 PM

I am working with the matrix-classes, and run into a problem. I need to multiply a SquareMatrix with a ColumnVector but I can't find an easy way to do so as multiplication isn't implemented on an IMatrix but only on the specific matrix-types themselves. The problem would be solved if we either had multiplication-methods that accept an IMatrix, or if all the matrix-types would inherit the Matrix type.

I'm also trying to get the Inverse() of a SymmetricMatrix, and this is currently not possible. There are several solutions to this problem I can think of:

  1. Implement the Inverse() method in SymmetricMatrix (and ISquareMatrix)
  2. Let SymmetricMatrix inherit SquareMatrix
  3. Implement Inverse() as an extension-method to ISquareMatrix

Thanks for the hard work and making this prommising library available for download!


May 11, 2009 at 7:38 AM

Hello Elmar, and thanks for writing. I think there are two issues at hand here: getting an inverse out of an ISymmetricMatrix and performing operations on an ISymmetricMatrix.

1. The lack of SymmetricMatrix.Inverse and a corresponding ISymmetricMatrix.Inverse is simply a serious bug that unfortunately got introduced in the clean-up before release. (We had such a method, but its implementation didn't take advantage of the symmetry. It was removed with the intent to replace it with a better algorithm, but when we decided to release before adding such an algorithm, we neglected to re-expose the old, less-than-optimally-efficient-but-still-worked-fine method.) I have filed a work item on it, and a fix should be available within a few days.

2. The inability to perform operations on interfaces is less easy to fix and has a long story behind it. Finding a good OO representation for matrices is ridiculously vexing. The matrix object model in Meta.Numerics went though at least four different incarnations trying to account for this. It did eventually become clear that using the class hierarchy to represent relationships between the various matrix classes wouldn't work. Here are some of the reasons:

a. There are multiple inheritance relationships, e.g. TridiagonalSymmetricMatrix should inherit from both TridiagonalMatrix and SymmetricMatrix. But the CLR doesn't allow multiple class inheritance. Multiple interface inheritance is allowed, so it works to use interfaces to represent relationships.

b. Classes carry not just a contracts but also an implementations. Thus, if SymmetricMatrix inherited from SquareMatrix, it would have whatever private storage SquareMatrix had, which isn't what we want. (It is possible to get around this by having an internal abstract layer around storage, but it get fairly awkward.) Interfaces carry only contracts, so this isn't a problem.

c. The override mechanism does not support changing the return type. For example, clearly SquareMatrix.Transpose() should return a SquareMatrix and SymmetricMatrix.Transpose() should retrun a SymmetricMatrix. But that means SymmetricMatrix.Transpose() can't be a virtual override of SquareMatrix.Transpose(). SymmetricMatrix.Transpose() could hide SquareMatrix.Transpose() using the "new" keyword, but then it wouldn't be called via late-bound binding when a SymmetricMatrix is passed to a routine that can handle any SquareMatrix. Interfaces get around this problem via explicit interface implementations, because they allow SymmetricMatrix to implement ISquareMatrix.Transpose() to return an ISquareMatrix seperately from ISymmetricMatrix.Transpose() to return an ISymmetricMatrix.

For all these reasons, we chose to use interfaces to repesent the hierarchy and classes to represent the implementations. This did have one very bad consequence which we haven't find a clean way around: since one cannot define operators on interfaces, there is no way to define an operation IMatrix * IMatrix, or any other operator overload involving pure interface types. (I have talked to the CLR people several times about this. It would have other good effects, too, for example allowing you constrain a template type to support operations like + or *. But it is not in the cards for the near future.) I see basically two options: (1) Live with it. The interfaces simply don't support operator, and users must code them explicitly. This is the current situation. (2) Define interface methods to replace the operators, e.g. IMatrix.Multiply(IMatrix). Coding matrix operations would be ugly because you couldn't use the intuitive operator notation, but it would still be possible without explictly specifying the algorithm. Do let me know which you prefer, or if you see a flaw in this analysis.

Cheers and thanks,

May 11, 2009 at 2:12 PM
Edited May 11, 2009 at 2:13 PM

Hello David,

Thanks for the quick reply and thank you for sharing the reasons behind the current design. I think of two ways to overcome the described problems (although I didn't try it out, so i'm potentially missing something):

1. Although it is not possible to implement operator overloads on interfaces, I think it is possible to implement overloads on operations between class-types and interfaces. This means we could implement an overload in the SquareMatrix-class that looks like this:

 public static IMatrix operator * (SquareMatrix M1, IMatrix M2)
          // Implementation

By giving every IMatrix-implementation an operator overload like this, we could overcome the current operator overload problem.

2. The other way around this problem would be to make a single "supermatrix" class (abstract class MatrixBase), that is inherited by all other Matrix-classes and only implements the operator-overloads between Matrix-objects (without accessing any internal storage). This wouldn't prohibit you from implementing a more efficient SquareMatrix * SquareMatrix - operator on the SquareMatrix-class, but would supply you with a basic 'fallback' operator-overload for any between-matrix operation.

I hope this reply is somewhat helpful. Thanks again!

Cheers, Elmar

May 11, 2009 at 7:01 PM

I think I might be misunderstanding what scenario you are trying to enable. If the scenario involves only interface types, e.g. you want to code

public static IMatrix MultiplyThreeMatrices (IMatrix A, IMatrix B, IMatrix C) {
return (A * B * C);

then your two suggestions won't help, because the compiler won't look for static operator methods on the SquareMatrix or any other class types.

If the scenario involves only class types, e.g. you want to code

SquareMatrix A = new SquareMatrix(3);
SymmetricMatrix B = new SymmetricMatrix(3);
// fill in A & B
SquareMatrix C = A * B;

then this should already be working for you. Is it not?

If the scenario is mixed, e.g. involves multiplying a SquareMatrix by an IMatrix, then your first suggestion would help.

This mixed scenario is not one we have fleshed out. We imagined that people would either be working in "class space" or in "interface space". "Class space" would be for people who are operating on matrices that their own code generates, so the exact type is known to the coder. "Interface space" would be for people who are operating on matrices that other code passes to theirs, so they would want to be able to get an inverse or eigenvalues without having to know the underlying type. (In this conext, you could say that our limitation is that operator notation is only supported in "class space".) If you have an example where you are driven to use a mixed scenario, I'd appreciate it if you could describe the requirements that led you there.

Thanks again!

May 11, 2009 at 8:10 PM

What I'm trying to achieve is exactly the "class types" scenario you describe. When I use your exact code, I'm getting an "Operator '*' cannot be applied to operands of type 'Meta.Numerics.Matrices.SquareMatrix' and 'Meta.Numerics.Matrices.SymmetricMatrix'"-error, because no multiplication overload is defined on this defined on these classtypes.

BTW: I just found out that the same scenario, involving a SymmetricalMatrix and a ColumnVector does work, because a multiplication-overload for IMatrix is defined on the ColumnVector-class.

Thanks for the help!


May 13, 2009 at 7:16 AM

You're absolutely right, Elmar. The latest check-in fixes both these issues. It exposes a SymmetricMatrix.Inverse method and defines overloads allowing +, -, and * operations between different matrix types.

I'm not thrilled with either of these fixes. The symmetric matrix inversion algorithm is not the LDL factorization algorithm that takes maximum advantage of the symmetry. The operations are done by defining explicit operator overloads for each pair of matrix types, a design that will lead to an N^2 explosion of operators as additional matrix types are added. (If you try to achive the same effect via interfaces, you will find that the compiler will complain that it can't decide which overload to use.) These are quick fixes that will enable functionality, but expect better fixes in the future.

Thanks for evaluating the product and helping us improve it! I look forward to continuing to interact with you in the future.