# Automated optimization of 0-1 matrix vector multiplication

Question:

Is there established procedure or theory for generating code that efficiently applies a matrix-vector multiplication, when the matrix is dense and filled with only zeros and ones? Ideally, the optimized code would make systematic use of previously computed information to reduce duplicated work.

In other words, I have a matrix $M$ and I want to do some precomputation based upon $M$, that will make computing $Mv$ as efficient as possible when I later receive the vector $v$.

$M$ is a rectangular dense binary matrix that is known at "compile time", whereas $v$ is an unknown real vector that is only known at "run time".

Example 1: (sliding window)

Let me use an easy small example to illustrate my point. Consider the matrix, $$M = \begin{bmatrix}1 & 1 & 1 & 1 & 1\\ & 1 & 1 & 1 & 1 & 1 \\ & & 1 & 1 & 1 & 1 & 1\\ & & & 1 & 1 & 1 & 1 & 1\end{bmatrix}.$$ Supposing we apply this matrix to a vector $v$ to get $w = Mv$. Then the entries of the result are, \begin{align} w_1 &= v_1 + v_2 + v_3 + v_4 + v_5\\ w_2 &= v_2 + v_3 + v_4 + v_5 + v_6\\ w_3 &= v_3 + v_4 + v_5 + v_6 + v_7\\ w_4 &= v_4 + v_5 + v_6 + v_7 + v_8 \end{align}

Doing a standard matrix-vector multiplication will compute exactly this way. However, a lot of this work is redundant. We could do the same matrix computation at less cost by keeping track of a "running total", and adding/subtracting to get the next number: \begin{align} w_1 &= v_1 + v_2 + v_3 + v_4 + v_5\\ w_2 &= w_1 + v_6 - v_1\\ w_3 &= w_2 + v_7 - v_2\\ w_4 &= w_3 + v_8 - v_3 \end{align}

Example 2: (hierarchical structure)

In the previous example, we could just keep track of a running total. However, usually one would need to create and store a tree of intermediate results. For example, consider $$M = \begin{bmatrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1 \\ & & & & 1 & 1 & 1 & 1\\ 1 & 1 \\ & & & & 1 & 1 \\ & & 1 & 1 \\ & & & & & & 1 & 1\end{bmatrix}$$ One could compute $w = Mv$ efficiently using a tree of intermediate results:

1. Compute $w_5$ and $w_7$, and add them to get $w_3$.
2. Compute $w_4$ and $w_6$, and add them to get $w_2$.
3. Add $w_2$ and $w_3$ to get $w_1$

The structure in the examples above is easy to see, but for the actual matrices I'm interested in, the structure is not so simple.

Example 3: (low rank)

To clear up some confusion, the matrices are generally not sparse. Specifically, a method solving this problem needs to be able to find efficient methods to apply matrices where large blocks are filled with ones. For example, consider

$$M = \begin{bmatrix}1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & & \\ 1 & 1 & 1 & 1 & & \\ 1 & 1 & 1 & 1 & & \end{bmatrix}.$$

This matrix can be decomposed as a difference of two rank-1 matrices, $$M = \begin{bmatrix}1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1\end{bmatrix} - \begin{bmatrix}& & & & & \\ & & & & & \\ & & & & 1 & 1 \\ & & & & 1 & 1 \\ & & & & 1 & 1\end{bmatrix}$$

so its action on a vector $w := Mv$ can be computed efficiently by, \begin{align} w_1 &= v_1 + v_2 + v_3 + v_4 + v_5 + v_6 \\ w_2 &= w_1 \\ w_3 &= w_2 - v_5 - v_6 \\ w_4 &= w_3 \\ w_5 &= w_4. \end{align}

Motivation:

I'm working on a numerical method for some image processing, and there are several large dense $0-1$ matrices with different structures that are fixed for all time. Later these matrices will need to be applied to many unknown vectors $v_i$ that will depend on the user's input. Right now I'm using pencil-and-paper to come up with efficient code on a per-matrix basis, but I'm wondering if the process can be automated.

Edit: (postscript)

All of the answers here so far (as of 9/5/15) are interesting, but none answer the question as satisfactorily as I had hoped. Probably it turns out that this is a hard research question, and no one knows a good answer.

Since time has run out I am awarding the bounty to EvilJS's answer since it addresses the right question. However, I wish the answer contained more clear and detailed explanations.

tranisstor's answer makes a connection between this question and the Online Boolean Matrix-Vector Multiplication (OMv) problem, but the connection is not exactly what this question is asking. In particular, the following assumption doesn't really fit (bold emphasis mine),

Now assume that for all $n \leq n_0$ and all $n \times n$ matrices $M$ we know an algorithm $A_{n,M}$, that for all vectors $v$ computes $Mv$ in truly subquadratic time, i.e. in time $O(n^{2 - \varepsilon})$ for some $\varepsilon > 0$.

Whether or not there exist subquadratic algorithms for all matrices is orthogonal to the question of finding an algorithm for a specific matrix that is as fast as possible. Most 0-1 matrices look like random noise and (if I were to guess) probably don't have subquadratic algorithms. However, the fact that there are really bad matrices out there doesn't prevent me from finding a fast algorithm on a good matrix, for example, a "sliding window" matrix.

vzn's answers, first answer, second answer are interesting (and in my opinion don't deserve so many downvotes), but they don't apply to the question for reasons discussed in the comments there.

• If your matrix is of this form, TDMA this is band matrix, Thomas algorithm. Not yet 0-1 but this feature should be exploited. – Evil Aug 26 '15 at 19:08
• @EvilJS the matrix just happens to be banded for the particular example. In general it will not be banded. I have added another example that is not banded. – Nick Alger Aug 26 '15 at 23:03
• You have a lot of constant matrices N x M which are binary, real vectors, and want to precompute optimal execution path during preprocessing stage per instance? The output of such operation is code with hardcoded operations per matrix and you want method to do so? By per instance I mean per matrix. Just checking. – Evil Sep 2 '15 at 16:59
• @EvilJS This question is about the situation where there is one known binary matrix $M$, which will be applied to many unknown real vectors $v_i$ at a later time. Based on $M$ only, we wish to precompute a code that will apply $M$ as efficiently as possible, so that later when we receive the $v_i$, we can compute $M v_i$ as fast as possible. In the particular application that motivates this question I have a handful of binary matrices like this (12 actually) that are fixed for all time whereas the vectors $v_i$ are unpredictable and depend on input from the user of the program. – Nick Alger Sep 2 '15 at 17:31
• Over the field of two elements, the problem of computing the minimum XOR-gate circuit that simulates a given linear transformation is NP-hard. See cstheory.stackexchange.com/a/32272/225 – Ryan Williams Sep 26 '15 at 1:13

If it is possible try to exploit banded tridiagonal nature of matrix.
Otherwise if the matrix contains only a constant number of distinct values (which surely is being binary), you should try Mailman algorithm (by Edo Liberty, Steven W. Zucker In Yale university technical report #1402): optimized over finite dictionary
Common Subexpression Elimination is know for some time like Multiple Constant Multiplication, but going down to gate-level is an option - patterns used here could be used separately as solution or merged with other methods, the paper for this "Improving Common Sub-expression Elimination Algorithm with A New Gate-Level Delay Computing Method" by Ning Wu, Xiaoqiang Zhang, Yunfei Ye, and Lidong Lan published in "Proceedings of the World Congress on Engineering and Computer Science 2013 Vol II WCECS 2013, 23-25 October, 2013, San Francisco, USA" Gate level CSE

There is also crude but working method, to generate symbolic matrix with constants, vector with variables and plug it to Static Single Assingment (SSA) from compilers, which automates process of writing matrices by hand.

new algorithm prototype
What you did with running sum: \begin{align} w_1 &= v_1 + v_2 + v_3 + v_4 + v_5 \\ w_2 &= w_1 + v_6 - v_1 \\ w_3 &= w_2 + v_7 - v_2 \\ w_4 &= w_3 + v_8 - v_3 \end{align}
Gives 10 operations, and with my initial idea to use Thomas it is equivalent.
For now I am still writing and testing new algorithm, also runtimes are nasty, but first test result gave me surprising answer:

\begin{align} tmp_1 &= v_2 + v_3 + v_4 + v_5 \\ w_1 &= v_1 + tmp_1 \\ w_2 &= tmp_1 + v_6 \\ w_3 &= w_2 + v_7 - v_2 \\ w_4 &= w_3 + v_8 - v_3 \end{align}
Which gives 9 operations, defining them as + or - is 1, and = is 0.

\begin{align} w_1 &= v_1 + v_2 + v_3 + v_4 + v_5 + v_6 \\ w_2 &= w_1 \\ w_3 &= w_2 - v_5 - v_6 \\ w_4 &= w_3 \\ w_5 &= w_4. \end{align}

This gives 7 operations, my algorithm result gave:
\begin{align} tmp_1 &= v_1 + v_2 + v_3 + v_4 \\ tmp_2 &= v_5 + v_6 \\ w_1 &= tmp_1 + tmp_2 \\ w_2 &= w_1 \\ w_3 &= w_2 - tmp_2 \\ w_4 &= w_3 \\ w_5 &= w_4. \end{align}
Which gives 6 operations For now I can tell that I am using Hamming distance, & and | bitwise operations, counting usages and making something like Cocke–Younger–Kasami (CYK) - "a parsing algorithm for context-free grammars, named after its inventors, John Cocke, Daniel Younger and Tadao Kasami. It employs bottom-up parsing and dynamic programming." - from Wikipedia This is the same technique I use to build blocks of variables.

• (re rev5) plz give ref on "evergreen method". also, what is SSA? CYK dynamic algorithm? – vzn Sep 5 '15 at 18:29
• I have awarded the bounty to this answer, and explained why in an edit to my original question. – Nick Alger Sep 6 '15 at 1:35

This is related to an open research question, which is known as the "Online Boolean Matrix-Vector Multiplication (OMv) problem". This problem reads as follows (see [1]): Given a binary $n \times n$ matrix $M$ and $n$ binary column vectors $v_1, \dots, v_n$, we need to compute $M v_i$ before $v_{i+1}$ arrives.

Notice that the problem from the question is somewhat more general: It allows for $m \times n$ matrices and real-valued vectors. Observe that the problem with $n \times n$ matrices and Boolean vectors is "easier", as it presents a special case.

Clearly, the naïve algorithm for the Online Boolean Matrix-Vector Multiplication problem (which just uses standard matrix-vector-multipliction) takes time $O(n^3)$. There is a conjecture (see e.g. [1]) that this cannot be done truly faster than $O(n^3)$. (In more detail, this conjecture goes as follows: There exists no truly subcubic algorithm, which solves the Online Boolean Matrix-Vector Multiplication Problem, i.e. there is no algorithm with running time $O(n^{3 - \varepsilon})$ for $\varepsilon > 0$).

It is known that Williams's algorithm solves this problem in time $O(n^3 / \log^2 n)$. See [2] for more details.

It would be a breakthrough in the area of conditional lower bounds, if one could prove or disprove the above conjecture.

[1] Unifying and Strengthening Hardness for Dynamic Problems via an Online Matrix-Vector Multiplication Conjecture. by Henzinger, Krinninger, Nanongkai and Saranurak
[ http://eprints.cs.univie.ac.at/4351/1/OMv_conjecture.pdf ]

[2] Matrix-vector multiplication in sub-quadratic time: (some preprocessing required). by Williams
[ http://dl.acm.org/citation.cfm?id=1283383.1283490 ]

Update

One of the questions in the comments was as follows: We know $M$ at compile time. Can't we adjust our algorithm to suit $M$, so the OMv problem (conjecture) does not apply? We will see that this is not the case, unless the OMv conjecture fails.

The proof idea is simple: Assume we could give fast algorithms for all matrices up to some certain size (e.g. distinguishing all possible cases). After this certain size we use divide and conquer.

Here are the details:
Fix some $n_0 \in \mathbb{N}$, which (without loss of generality) is a power of 2 and bigger than 2. Now assume that for all $n \leq n_0$ and all $n \times n$ matrices $M$ we know an algorithm $A_{n,M}$, that for all vectors $v$ computes $Mv$ in truly subquadratic time, i.e. in time $O(n^{2 - \varepsilon})$ for some $\varepsilon > 0$. (Notice that this allows an individual algorithm for each matrix up to size $n_0 \times n_0$.)

Now we will solve OMv in truly subcubic time:
Given a binary matrix $M$ of size $n \times n$, where $n = 2^k$ for some $k$ and $n > n_0$, we use a divide and conquer strategy. We divide $M$ into four submatrices $M_1, M_2, M_3, M_4$ of sizes $2^{k-1} \times 2^{k-1}$. If $2^{k-1} \leq n_0$, then we use algorithm $A_{2^{k-1},M_i}$, otherwise, we recurse. (As $n_0$ is some fixed number, we can pick the correct algorithm in constant time.)

Notice that we will need at most $O(\log n)$ recursion steps. Also, for $n$ vectors $v_1, \dots, v_n$, we will $n$ computations. Thus, to process all matrix-vector multiplications we will need a total computation time of $O(n^{3 - \varepsilon} \log n)$.

It is well known that the logarithm grows slower than any polynomial (in particular slower than any root). Fixing some $\tilde \varepsilon > 0$ with $\tilde \varepsilon < \varepsilon$, we see that our total computation is running in truly subcubic time (in particular, in time $O(n^{3 - \tilde \varepsilon})$). Thus, the OMv conjecture would be wrong.

(If $M$ has size $m \times n$ and $m$ and $n$ are not powers of 2, then the bounds on the running times still apply, as we could just increase $n$ and $m$ to the next powers of 2.)

Conclusion: If you could make use of case distinctions on the input matrices to derive fast algorithms, then you could improve the OMv conjecture.

• As pointed out by author and vzn, this is not the case, vector is not binary, Matrix is not necessary N x N, and author wants to precalculate operations, and there is no need for online processing. Based on conjecture is not enough. Both papers are irrelevant to question. The case here is to precompute constant matrix to provide minimal number of operations. There will be possible different approaches for full, banded, symmetric cases. – Evil Sep 2 '15 at 16:55
• @EvilJS: If you allow any M x N matrix and real-valued vectors, then the problem just gets harder than the one I gave in the answer (i.e. the Online Boolean Matrix-Vector Multiplication will be a special case). If you could solve the more general problem truly faster than O(n^3), then you would also make an improvement on the conjecture (which would be big news!). Furthermore, the author says in a comment to the question that the vectors are initially unknown. If you knew all vectors beforehand, you could just use fast matrix multiplication (e.g. a version of Strassen's algorithm). – tranisstor Sep 2 '15 at 18:12
• I just pointed authors case "real vector". Look at Thomas matrix - only special case of matrices in O(n). I do not imply general case. And if Matrix are constant and vectors are known you hardcode answer not implement Strassen ;( – Evil Sep 2 '15 at 18:31
• @EvilJS: I am not sure I completely understand what you are trying to say. Of course, for special types of matrices like the Thomas matrix you can get a significant speed up, but in general this is harder. Maybe I should also point out that the problem I introduced does consider a preprocessing step (before any vector arrives). If you could tell me how to systematically "hardcode" your algorithm for any matrix that I give you, you could also improve on the conjecture (since you could implement this hardcoding step as a preprocessing step of an algorithm). – tranisstor Sep 2 '15 at 18:43
• agreed it works; however the 2nd ref by williams does not seem to consider binary matrices at all in particular. fyi he has slides here – vzn Sep 2 '15 at 18:55

this is essentially research-level CS, the problem is studied in at least two guises, one of multiplication of sparse matrices (example paper just cited), and also the special case of "binary sparse matrices" is also studied. the 2nd case is known to be related to optimizing straight-line programs. minimal programs may also be like DAGs with two types of "gates", addition and multiplication, so some circuit minimization literature may connect with this, and possibly "off the shelf" software could be adapted for the purpose. here is a specific ref on the 2nd case and also the same question on cstheory with some basic initial empirical study.

• I only skimmed the links, but so far I'm skeptical for two reasons. 1) The matrices are not sparse. Often the number of nonzeros is around the same as the number of zeros, yet there is often a way to come up with $O(n)$ algorithms for matrices with $O(n^2)$ nonzeros, based on the patterns. 2) These papers seem to be based on only using addition and multiplication as the basic binary "gates", but in my experience it is important to also use subtraction. This would require including a subtraction binary gate, or a "multiply by -1" unary gate. – Nick Alger Aug 28 '15 at 21:55
• the refs are on, as the titles indicate, sparse matrices. maybe you have some different definition than those in the papers? if you are sensitive to an exact definition of sparsity (most are roughly correlated/ nearly interchangeable) it should be stated in the question. – vzn Aug 31 '15 at 2:15
• The matrices I'm interested in are dense matrices. By the way, even though I don't think this fully addresses my question, I appreciate the answer. – Nick Alger Aug 31 '15 at 4:38
• ok, sorry! got mixed up, didnt realize the exact question. on cursory look your example #2 has less than ½ fill & looked "sparse" to me & figured some of the sparse theory would be at least somewhat applicable. basically the more dense the matrix, the less the operation can be optimized, so probably most of the theory about this type of optimization is oriented around sparse matrices. – vzn Aug 31 '15 at 15:03

not sure if this problem has been studied exactly but this research is related & seems a reasonable lead/ start. it looks at hypergraph decomposition for sparse matrix multiplication. binary matrices are a special case of this approach. this approach will find more optimal strategies than the "straight" multiplication method. further optimizations (within this framework) might be possible based on the binary matrix property.

• I don't see what this has to do with the question. That paper is about partitioning matrix multiplication among a distributed system, for parallel computation, to minimize the amount of inter-processor communication. What does that have to do with this question? The question does not seem to mention anything about parallel computation or inter-processor communication. I encourage you to edit your answer to make the connection more explicit. – D.W. Aug 27 '15 at 19:01
• afaik its the same problem and minimizing the parallel computation also minimizes a single processor implementation of the same calculations. at least, the questioner did not rule out parallel implementations. – vzn Aug 27 '15 at 19:39
• Thank you for the link. However I am skeptical of the method for this problem since it does not take advantage of the fact that the matrix entries contain only zeros and ones, whereas this property is very important, as far as I can tell. For example, the "running total" algorithm in the first example will only work if all the nonzero entries in a given column of the matrix have the same value. – Nick Alger Aug 28 '15 at 4:49
• NA your observation/ objection is addressed in the answer. further optimization is probably possible using the 0/1 property; this method seems to minimize total # of addition/ multiplication operations under the guise of parallelization. the addition/ multiplication operations can also be seen as "gates" in a DAG and the technique is minimizing gates. the substantial complexity of the paper reveals some of the inherent deeper/ substantial complexity of this optimization process. as stated answer is not intended to be definitive on this difficult problem, just "better than nothing". – vzn Aug 28 '15 at 14:16