C++ Matrix class hierarchy

Question

Should a matrix software library have a root class (e.g., MatrixBase) from which more specialized (or more constrained) matrix classes (e.g., SparseMatrix, UpperTriangluarMatrix, etc.) derive? If so, should the derived classes be derived publicly/protectively/privately? If not, should they be composed with a implementation class encapsulating common functionality and be otherwise unrelated? Something else?

I was having a conversation about this with a software developer colleague (I am not per se) who mentioned that it is a common programming design mistake to derive a more restricted class from a more general one (e.g., he used the example of how it was not a good idea to derive a Circle class from an Ellipse class as similar to the matrix design issue) even when it is true that a SparseMatrix "IS A" MatrixBase. The interface presented by both the base and derived classes should be the same for basic operations; for specialized operations, a derived class would have additional functionality that might not be possible to implement for an arbitrary MatrixBase object. For example, we can compute the cholesky decomposition only for a PositiveDefiniteMatrix class object; however, multiplication by a scalar should work the same way for both the base and derived classes. Also, even if the underlying data storage implementation differs the operator()(int,int) should work as expected for any type of matrix class.

I have started looking at a few open-source matrix libraries and it appears like this is kind of a mixed bag (or maybe I'm looking at a mixed bag of libraries). I am planning on helping out with a refactoring of a math library where this has been a point of contention and I'd like to have opinions (that is unless there really is an objective right answer to this question) as to what design philosophy would be best and what are the pros and cons to any reasonable approach.

Answer 1

The problem with the Circle subclass of an Ellipse (or a Square subclass of a Rectangle) occurs when you can modify one dimension per the Ellipse interface, so that the circle is no longer a circle (and the square no longer square).

If you only allow nonmodifiable matrices then you are safe, and you can structure your type hierarchy in the natural way.

Answer 2

hehehe. At first I read that your friend was saying that a Circle should be an Ellipse and wrote a long tirade about why they were full of it.

You should listen to your friend, except that I hope they're not saying that a SparseMatrix "is-a" MatrixBase. The term means different things in the real world vs. the modeling world. In the modeling world, "is-a" means following the Liskov Substitution Principle (look it up!). Alternatively it means that SparseMatrix must follow the contract of MatrixBase in that member functions must not require any extra preconditions and must meet no less postconditions.

I don't know exactly how this applies to the matrix issue but if you look into the terms I used in the previous paragraph (LSP and Design by Contract) then you should be well on your way toward learning the answer to your problem.

One way that might apply in your case is to take the various commonalities across your hierarchy and make them abstract interfaces. Then inherit from these interfaces in those classes that respond to them correctly. This would allow you to write functions that should allow common use and yet still retain separation where there is too much variation.

Answer 3

This is a nice question, but I am not yet sure what the metrics are by which you want to evaluate this.

For what it is worth, the one Matrix library I currently use the most is Armadillo does have a common Base object using the curiously recurring remplate pattern. I believe Eigen (another recent and heavily templated Matrix library) does the same.

Answer 4

The main problem to watch out with inheritance-based designs like this is SLICING.

Let's say MatrixBase defines a non-virtual assignment operator. It copies all the data members common to all matrix subclasses. Your SparseMatrix class defines additional data members. Now what happens when we write this?

SparseMatrix sm(...);
MatrixBase& bm = sm;
bm = some_dense_matrix;

This code makes little sense (trying to assign a DenseMatrix to a SparseMatrix directly through an operator defined in the base class) and is prone to all kinds of nasty slicing behavior, yet is a vulnerable aspect of such code and there's a very big possibility of this happening somewhere down the line if you provide assignment operators accesible through MatrixBase*/MatrixBase&. Even when we have this:

SparseMatrix sm(...);
MatrixBase& bm = sm;
bm = some_other_sparase_matrix;

... we still have slicing issues due to the assignment operator being non-virtual. Without the inheritance of a common base class, we can provide assignment operators to meaningfully copy a dense matrix to a sparse one, but trying to do this through a common base class is prone to all kinds of issues.

Assignment operators should generally be avoided for base classes! Imagine a case where Dog and Cat inherit from Mammal and Mammal provided an assignment operator, virtual or not. It would imply that we can assign Dogs to Cats, which makes no sense and even if the operator were virtual, it would be difficult to provide any kind of meaningful behavior to assign mammals to other mammals.

Let's say we try to improve the situation by implementing the assignment operator in Dog so that it can only be assigned other dogs. Now what happens when we inherit from Dog to create Chihuahua and Doberman? We should not be able to assign Chihuahuas to Dobermans, and so the original case recursively repeats itself until you're sure you've arrived at the leaf nodes of an inheritance hierarchy (it's a shame C++ does not have a final keyword to prevent any further inheritance).

The same problem is apparent with the common Circle inherits Ellipse example. The circle might require width and height to match: that's an invariant the class wants to maintain, yet anyone can simply get a base pointer (Ellipse*) pointing to a Circle object and violate that rule.

If in doubt, avoid inheritance as that is one of the most misused features of C++ and any language which supports object-oriented programming in general. You can try to work around the problem by providing runtime mechanisms to determine the type of a subclass being assigned to another subclass and only allow for matching types, but now you're doing a whole lot of extra work and incurring runtime overhead. It's better to avoid assignment operators all together for inheritance hierarchies and rely on methods like clone to produce copies (prototype pattern).

Thus, if you choose to make an inheritance hierarchy out of your matrix classes, you should think carefully as to whether the (most likely short-term) advantages of inheriting outweigh the long-term disadvantages. You should also be sure to avoid all potential cases where slicing can occur, which might be very difficult to do for a matrix library without compromising its usability and efficiency.