# What's the proper class hierarchy to use for matrices?

I'm working with particular classes of matrices which admit fast operations without resorting to dense linear algebra. To keep the heierarchy simple, I'm currently doing the following (greatly simplified):

``````class Matrix:
def mvm() # matrix-vector multiplication
def is_square()
def is_psd()
def is_symm()
def eigenvalue_bound()

class Diag(Matrix):
# represents a diagonal matrix with main diagonal v
# requires O(n) memory and time for all Matrix operations
# where n is the size of v
def __init__(self, v)

class Circulant(Matrix):
# represents a circulant matrix with top row v
# requires O(n lg n) memory and time for all Matrix operations
# where n is the size of v
def __init__(self, v)

class Kronecker(Matrix):
# represents the Kronecker product of two matrices A and B
# requires O(memory(A) + memory(B)) memory for sparse storage
# instead of the dense version, of size O(memory(A) * memory(B))
# Similarly, a Matrix operation taking op(A) and op(B) time
# can be done in O(size(A) * op(B) + size(B) * op(A)).
def __init__(self, A, B)

class Composition(Matrix):
# represents the composition product of two matrices A and B, AB
# requires O(memory(A) + memory(B)) memory
# might save time if op(A) + op(B) < op(AB)
def __init__(self, A, B)
``````

As you can see, I care about a lot of properties of the matrices - positive semi-definiteness, symmetry, eigenvalue information. These properties are available only if the matrix meets certain mathematical properties, and they translate nontrivially between the matrices.

There are other matrices and properties I'm considering, but hopefully this gets the point across: the methods I'd like to call on an instance of each of these is heavily dependent on the context in which the matrices are created.

I've adopted the "simple" solution of having the base contain the union of all properties I might ever want, and then I perform various runtime sanity checks if a method is called in a certain context.

At the other end of the spectrum, this might be a solution for multiple inheritance and mixins, but this would result in a combinatorially large number of subclasses to handle! What's the best approach here?

For clarification: these matrices are immutable. For consistency, let's stick to python. Another important feature is that under different constraints, the in-memory representation differs between each subtype. In other words, I'm really just treating these matrices as linear operators.

• is the matrix immutable? if not then the matrix may change from identity to fully arbitrary. – ratchet freak Jan 30 '17 at 15:51
• I'm not crazy about inheritance for this. Would suggest you favor composition. en.wikipedia.org/wiki/Composition_over_inheritance – Brad Thomas Jan 30 '17 at 16:01
• @gardenhead sure it does. The immutability is just a contract. Python just doesn't hold you to it (this is true for lots of other language features of python). Either way, it's beside the point of the question. – VF1 Jan 30 '17 at 19:20
• @VF1 Well, to get to the point of the question, you shouldn't be using inheritance at all. Are you coming from a math background? – gardenhead Jan 30 '17 at 20:24
• If you want an example of code following a highly typed strategy, you can look at the linalg package for Julia. – Andrew Piliser Jan 31 '17 at 18:28

Given your matrixes are immutable, it is indeed possible to have a general base class `Matrix` which defines all the operations on matrices you need, and a few specific subtypes of this class as a form of optimization, when you know you need this. However, before doing so, I would really think twice if the kind of additional complexity you add to your program is really worth the hassle, and if you are not creating a maintenance nightmare by doing so.

For example, a subtype `Diag` derived from `Matrix` can eventually make sense, since it can be implemented with much better memory & time behaviour for some of its typical operations than arbitrary matrices. I guess it can bring benefits if you are having use cases where you need lots of them.

The need for a subtype `SquareMatrix` however is debatable, since it won't bring you any significant memory savings. Maybe it is true what you wrote in "various implementations can be greatly simplified if the matrices are square", but are those features or functions you really need? For example, matrix multiplication does not fall into this category, the implementation effort is essentially the same for square and non-square matrices.

Similar considerations should be done for subtypes like "Circulant" or "Kronecker" or "Composition". Having a specific subtype for these lets you implement some operations like matrix inversion or test for PSD in an optimized way. Note, however, that other operations like standard matrix addition or extracting single entries will loose efficiency. So I would only provide such specific subtypes if you have a real use case for them which needs to be optimized and where you are sure it is worth the hassle.

So instead of thinking about how to build as many subtypes as possible, I recommend you ask yourself how to use as few subtypes as possible, and optimize only when you really run into measurable performance or memory problems. Optimization is not an end in itself, it is a means to an end, and if you do not have a good reason why you need these optimizations, you should consider not to optimize at all.

• Ok, this seems to be the most reasonable advice - keep with the base containing a union of operations, then use subtypes judiciously. See my reply on my question to your comment regarding recursively defined matrices. – VF1 Jan 30 '17 at 23:23
• @VF1: your quiet change of "Dense" to "Hankel" invalidated parts of my previous answer. But FWIW, I tried to change my answer to what you may have in mind, please leave a note when I got you wrong. – Doc Brown Jan 31 '17 at 6:51
• My apologies, I thought changing to Hankel would simplify the number of cases you had to look at (which it seems to have done). Also, I wanted to give an example of a structured yet non-sparse matrix. I think your answer captures the architectural question here, namely that there's a tradeoff between complexity and speed that we can tune by using the two-level Matrix hierarchy and a well-chosen number of subtypes. I'll keep the question open for a bit longer, though, in case someone has a deep insight. – VF1 Jan 31 '17 at 13:58

Using a class hierarchy is not generally advisable in such cases, since your subtypes encode specific data constraints.

If your matrices are mutable (i.e. elements of the matrix can be updated), then you run into the circle–ellipse problem. E.g. given a diagonal matrix, I could change an element so that the object no longer represents a diagonal matrix. However, once constructed, an object cannot change its type.

The next problem is that your constraints may not be exclusive. This depends on what exact constraints you are trying to encode. Unless you are using a mixing system, this leads to a combinatorial explosion of types. If you are using a flexible mixin system, you still can't select the correct type for a matrix in your program without knowing the properties of the runtime data.

As a consequence, the only reasonable solution is to have a single `Matrix` type without any particular subtypes.

How can we then deal with operations that only make sense when the matrix fulfils a particular constraint? If these are part of the general `Matrix` type, they may cause runtime exceptions. This can be made more obvious by providing methods that check the constraint. If the constraint is met, they return an adapter object providing the restricted operation. If the constraint is not met, they return null. In C++, this would be used as:

``````// instead of: matrix.try_operation_only_for_diagonals()

if (auto matrix_as_diagonal = matrix.as_diagonal()) {
matrix_as_diagonal->operation_only_for_diagonals();
}
``````

In Python:

``````matrix_as_diagonal = matrix.as_diagonal()
if matrix_as_diagonal is not None:
matrix_as_diagonal.operation_only_for_diagonals()
``````

The object representing this constraint has only a single member field: the original matrix for which it supplies additional behaviour.

If your matrices are immutable or change rarely, you can cache the results of the constraint checks in the `Matrix` object. Whether all constraints should be checked during matrix construction or lazily depends on your application. If you can statically prove that a matrix will satisfy certain constraints, allowing these to be set in the constructor would be best. How that can be done elegantly (i.e. without a combinatorial explosion of constructor signatures) depends on your programming language.

How can we deal with operations that are available for all matrices but are implemented differently when a certain constraint is met? By the same mechanism. We need to test for the various constraints at runtime. But it would be annoying if the constraints would be checked again and again for each execution of the operation. In such cases, you can use the strategy pattern. The first time you execute such an operation, you select an appropriate execution strategy, and can later use it directly when that operation is executed again on the same matrix.

• If it helps clarify/focus your answer, you can assume immutable matrices + Python. I'll edit the question accordingly. I'm still not fully understanding what it means to have no subtypes. Differenly-constrained matrices are constructed differently and have different representations in memory. The whole reason I'm not just using `numpy` is because, for instance, I can do all operations I want with only a vector representing the main diagonal of the diagonal matrix. The same argument follows for differently-structured matrices. – VF1 Jan 30 '17 at 18:38
• +1 for encoding the constraint as a type when using immutable objects. It's a very elegant pattern and frees you from having to re-check the constraint over and over again or tracking a separate boolean. – Doval Jan 31 '17 at 1:17
• @Doval am I misunderstanding something? this answer recommends getting rid of the types. – VF1 Jan 31 '17 at 13:56

My answer comes late, but I believe it is still important to share my insights with you. Technically it is not my insight, but instead the insights learned from OpenCV.

TL; DR:

• Doc Brown's answer favors a no-go decision by default, and urges you to think critically whether the potential benefits are exploitable in practice or not.
• My answer explains how existing libraries handle the same decision. But my viewpoint is really the same: don't go.

These other matrices are best expressed (in the programming language's type system) as Matrix Expressions. In other words, they are special objects that can be evaluated to give a matrix. However, each special object also has queryable properties, which can be used by linear operator implementations to perform the task, without having to inspect every `(i, j)` of the evaluated matrix.

Matrix expression is a generalization of a matrix. In programming terms, a matrix expression is a discriminated union (sum type) of:

• Matrix expression
• Fully-populated matrices
• Any one of the special matrices - a list that is potentially extensible
• Unary operator applied to a (child) matrix expression
• Binary operator applied to two (left, right) matrix expressions

The following is the motivation of why expression evaluation on special matrices shall be deferred.

• The programmer-user assigns the product of A, B to C
`C = A * B`
• The programmer-user assigns the product of C, D to E
`E = C * D`
• Matrix `A` is typical (dense, not having special properties)
• Matrices `B` and `D` are special.
• At the time of assignment of `(A * B)` to `C`, because one of the argument ('A') is typical, it appears as if a full evaluation of the matrix product is needed.
• But, once it is known that `C` is to be multiplied with `D`, which is special, it has already missed the chance to carry out an optimization, which is `A * (B * D)`.

Therefore, if the evaluation of matrix expressions are not deferred, a lot of optimization opportunities would have been forgone due to the architectural decision.

While matrix multiplication is associative, there might be other operators that aren't. That has to be kept in mind when implementing matrix expressions.

Familiarity with Einstein summation would be helpful in a discussion: http://mathworld.wolfram.com/MatrixMultiplication.html

Matrix properties (such as symmetry) are entirely separate. Properties behave like predicates. A matrix expression can be passed into a predicate to get a `true/false` result.

However, for the performance-conscious users, each predicate will have several performance characteristics:

• Constant case: `is_symm ( MatExpr )` can be returned in `O(1)` time
• Linear case: `is_symm ( MatExpr )` can be returned in `O(M+N)` time
• Fully evaluated case: `is_symm ( MatExpr )` can be returned in `O(MN)` time
• Cacheability: the given `MatExpr` is found to be immutable, and fortunately its symmetry has already been computed or cached somehow, so that `is_symm ( MatExpr )` can return the result in `O(1)` time.

Algorithms that need to decide whether to make use of the symmetry property will have to query the performance characteristics first, and then decide the evaluation strategy.

Sounds complicated? Yes, so don't go. It's painful.

This point of view has been implemented in various matrix programming packages, so I am not bringing up anything new. Also, implementation is presumably in the hands of Profs and PhDs. I am not qualified to say anything as to how they are actually designed and implemented in software libraries.

• Great answer, thanks for the insights. Can you link some of the libraries you refer to? Also, I suppose in my case my rationale was that the user decides what should be kept evaluated or not. In particular, if you have `SumMatrix(Composition(A,B), Toeplitz(C)`, you're telling the internal representation to be that entire "deferred" tree. The job of deciding what should be simplified is more of a CAS problem, and in my case the composability achieves what I want. – VF1 Feb 3 '17 at 23:21
• @VF1 I have to break the bad news: OpenCV does not actually implement the full glory. It tried, but it's too hard. (It has `MatExpr` but it still evaluates matrix operations sequentially, and no support for special matrices.) As I said, I don't know the details of the other libraries. – rwong Feb 3 '17 at 23:24