An Efficient Matrix Multiplication Algorithm for Accelerating Inference in Binary and Ternary Neural Networks

In summary, our contributions are as follows:

• •

  (Section 3) We provide a formal definition of the vector-ternary-matrix multiplication problem and demonstrate its reduction to the vector-binary-matrix multiplication problem.

• •

  (Section 4) Following the observation that weight metrics are fixed once a model is trained, we preprocess them and construct indices that enable the development of our efficient algorithms while reducing the memory requirements by a logarithmic factor.

• •

  (Section 5) We introduce the RSR algorithm with a time complexity of $O\left(\frac{n^{2}}{\log(n)-\log(\log(n))}\right)$ for vector-binary-matrix multiplication for $n$ by $n$ matrices. We further optimize this algorithm and introduce RSR++ that achieves a faster running time of $O(\frac{n^{2}}{\log(n)})$.

• •

  (Sections 6 and 7) We conduct various experiments to demonstrate the applicability of our algorithms for matrix multiplication using different implementation configurations. We show that we can achieve up to 29x faster run time and 6x less space usage on matrix multiplication task. Finally, we discuss the advantages and limitations of our algorithm in Section 7.

We denote vectors by $\vec{v}$. We denote matrices by capital letters. For example, $A\in E^{n\times m}$ is a matrix of size $n\times m$ with elements belonging to a set of permissible values, $E$. Specifically, if $E=\{0,1\}$ (resp. $E=\{-1,0,1\}$), then $A\in E^{n\times m}$ is denoted as a binary (resp. ternary) matrix.

Let $\Sigma_{n}$ denote the set of all bijective permutations functions $\sigma:\{1,2,\dots,n\}\rightarrow\{1,2,\dots,n\}$. For a vector $\vec{v}$, the permuted vector under $\sigma$ is represented as $\pi_{\sigma}(\vec{v})$, where each element is repositioned such that $\pi_{\sigma}(\vec{v})[i]=\vec{v}[\sigma(i)]$, moving the element $\sigma(i)$ in $\vec{v}$ to position $i$ of the permuted vector. We extend $\pi_{\sigma}$ to matrices by permuting their rows. When $\sigma$ is clear by the context, we may simplify $\pi_{\sigma}$ to $\pi$. Let $\mathcal{L}^{\mathbb{N}}_{n}$ denote the set of all ordered lists of length $n$ with entries in $\mathbb{N}$, and let $\mathcal{L}^{\mathbb{N}}_{\leq n}$ represent the set of sorted lists with length at most $n$. We use $A[:,j]$ to indicate the $j$-th column of a matrix $A$. Similarly, we use $A[i,:]$ to refer to the $i$-th row the matrix $A$ and $A[i,j]$ to indicate an element.

Given a vector $\vec{v}\in\mathbb{R}^{n}$ and a ternary matrix $A\in\{-1,0,1\}^{n\times n}$, our objective is to efficiently compute the product $\vec{v}\cdot A$. In this setting, the vector $\vec{v}$ serves as the activation output from the previous layer of the neural network, and $A$ is the weight matrix (See Figures 1 and 2).

Here, we focus on a single multiplication instance and aim to accelerate this operation beyond the standard naive vector-matrix multiplication complexity of $O(n^{2})$. We initially focus on square matrices of size $n\times n$ for clarity in exposition. However, our approach is applicable more broadly and extends to general matrices of dimensions $n\times m$ as discussed in Section  7.

A key observation is that the weight matrices are fixed after a model is trained, while the vector $\vec{v}$ is provided at inference time. It allows us to preprocess each weight matrix $A$ to enable faster vector-matrix product – hence more efficient inference.

Initially, we use the following lemma to reduce our problem into a Vector-Binary-Matrix product.

As a result, we focus on solving the following problem.

As a result of Lemma  3.1, any time guarantees established for this reduced version will apply equivalently to the original problem by a constant factor.

We introduce the Redundant Segment Reduction (RSR) algorithm, which optimizes matrix-vector multiplication by identifying and reusing redundant computations. The RSR algorithm begins by partitioning the matrix into smaller blocks, each with a row size of $n$ but a reduced column size.

Within each block, similar segments across different columns are identified. We will calculate the dot product for each unique segment only once. This approach reduces the time for computing the product of the vector and each matrix block. To maximize the length of similar segments across columns, we also apply a permutation to the matrix blocks. In the final step, we aggregate the results of the vector’s product with each matrix block to obtain the overall product. An overview of this approach is illustrated in Figure 3.

In Section § 4, we present the preprocessing phase of our algorithm, detailing the construction of essential data structures utilized during inference time. Subsequently, in Section § 5, we describe the inference-time operations, illustrating how our algorithm accelerates the multiplication process, achieving a logarithmic improvement by reducing the $O(n^{2})$ time complexity of the standard vector-matrix-multiplication to $O\left(\frac{n^{2}}{\log(n)}\right)$.

The first step of our preprocessing is Column Blocking – the process of partitioning consecutive columns of the original binary matrix $B$ to construct a set of smaller, compact matrices, each representing a distinct block of columns.

For example, consider the following binary matrix,

The set of $2$-Column Blocks of $B$ are as follows:

When $k$ is clear by the context, we use $B_{i}$ instead of $B^{[k]}_{i}$ to denote the $i$-th block.

The next step after column blocking is binary row ordering, where the rows of each column block are sorted according to a lexicographic order.

To further clarify the Binary Row Order, let us consider the following example, where $B_{i}$ is a block with two columns. The permutation function $\sigma_{B_{i}}=\langle 2,5,6,1,3,4\rangle$⁴⁴4$\sigma_{B_{i}}(j)$ is the $j$-th value in $\sigma_{B_{i}}$. For example, $\sigma_{B_{i}}(2)=5$. provides a Binary Row Order $\pi_{\sigma_{B_{i}}}(B_{i})$ of the matrix $B_{i}$.

The right expression results from $00_{2}<01_{2}<11_{2}$.

We define $Bin_{[k]}$ as a binary $2^{k}\times k$ matrix which has exactly one row for each binary value from $0$ to $2^{k}-1$, and it is also Binary Row Ordered. For example,

Next, we segment each binary row ordered matrix (the sorted column blocks) into groups of rows with the same binary value representation, enabling further matrix compaction.

Consider the Binary Row Ordered matrix $\pi_{\sigma_{B_{i}}}(B_{i})$ shown in Example 2. The Segmentation list of this matrix is $[1,4,6]$ because the first row is the beginning of $00$ rows. Then, $01$ starts from index $4$, and $11$ starts at index $6$.

For a Binary Row Ordered matrix of dimensions $m\times k$, we extend this concept to define Full Segmentation as the list of exactly $2^{k}$ elements, where the $j$th element of the list indicates the first index of a row that has the binary value $j$. For example, the Full Segmentation of matrix $\pi_{\sigma_{B_{i}}}(B_{i})$ in Example  2 is $[1,4,6,6]$. Since we don’t have any rows for $10$, we assign the same boundary as the next value. See Table 1 for an illustration.

From now on, we use the same notation $\mathcal{S}$ to denote the Full Segmentation for a Binary Row Ordered matrix.

Given the Full Segmentation of a Binary Row Ordered matrix, a compact representation of the matrix can be obtained by retaining only the first indices of each unique row value in the segmentation list. In the following sections, for simplicity, we use $L_{i}$ instead of $\mathcal{S}(B_{i})$ to show the Full Segmentation of the permuted Column Block $B_{i}$.

Putting all three steps together, the preprocessing algorithm is described in Algorithm  1. The algorithm starts by selecting a parameter $k$ such that $k\leq\log_{2}(n)$. Next, it computes the $k$-Column Blocks of $B$. Hereafter, we denote these blocks as $B_{i}$ rather than $B^{[k]}_{i}$ for $i\leq\lceil\frac{n}{k}\rceil$. For each $B_{i}$, the algorithm proceeds to calculate both the Binary Row Orders $\sigma_{B_{i}}$ and the Full Segmentations $L_{i}$.

Figure 3 (Preprocessing – the left figure) provides a visual example of the preprocessing steps.

To simplify the complexities, in the following, we provide the time and space complexities for an $n$ by $n$ weight matrix. Note that storing the weights matrix has a space complexity of $O(n^{2})$.

Recall that during the preprocessing time, for each column block $B_{i}$, we construct a permutation function $\sigma_{B_{i}}$ and a full segmentation list $L_{i}$ as its index.

At the inference time, given a vector $\vec{v}$, our objective is to efficiently compute the multiplication of $\vec{v}$ by each column block $B_{i}$, using $(\sigma_{B_{i}},L_{i})$.

Our first step is computing the Segmented Sums over the permuted vector $\pi_{\sigma_{B_{i}}}(\vec{v})$ – i.e., the sum for each interval in the permuted vector that corresponds to groups of similar rows in the binary matrix $B_{i}$.

For example, for the Full Segmentation defined for Example  2 and a given vector $\vec{v}=[3,2,4,5,9,1]$, we have:

Where $9=3+2+4,14=5+9,1=1$, and the third element is 0 because there is no $10_{2}$ in the segmentation.

Computing the Segmented Sums using Equation 3 requires first computing the permuted vector $\vec{v}_{\pi}$. Instead, to efficiently compute the Segmented Sums in place (without computing the permuted vector), we use Equation 5:

We are now ready to describe our algorithm RSR for computing $\vec{v}.B$ at the inference time. Given the input vector $\vec{v}$, for each column block $B_{i}$, in Step 1 we use the permutation function $\sigma_{B_{i}}$ and the full segmentation list $L_{i}$ to compute the Segmented Sum $SS_{L_{i},\sigma_{B_{i}}}(\vec{v})$ using Equation 5.

Then, we use Lemma 5.2 to calculate $\vec{v}\cdot B_{i}$.

Following Lemma 5.2, as the final step in the inference time, in Step 2 we calculate the product $SS_{L_{i},\sigma_{B_{i}}}(\vec{v})\cdot Bin_{[k]}$ for each $k$-Column Block $B_{i}$. Algorithm  2 shows a pseudo-code of this algorithm.

Having presented RSR, we now present our faster algorithm RSR++ that focuses on improving Step 2 of RSR. This step involves computing the product of a vector $\vec{u}$ of size $2^{k}$ and the binary matrix $Bin_{[k]}$ of size $2^{k}\times k$ (Line 5 of Algorithm  2).

The standard vector-matrix multiplication, in this case, leads to $O(k\cdot 2^{k})$ time. However, the special structure of matrix $Bin_{[k]}$ allows us to reduce it to $O(2^{k})$. Algorithm  3 shows the pseudo-code for this multiplication. See Figure 4 for a visualization of this approach.

Let $\vec{r}$ denote the final result of this product, $\vec{u}\cdot Bin_{[k]}$. We build $\vec{r}$ starting from the $k$-th element to the first one. (i) The $k$-th element can be calculated by summing up the elements in $\vec{u}$ with odd indices. For example, in RSR++ visualization in Figure 4, the sum of odd-index elements (orange cells) in the top row is 23. Then, (ii) we calculate a vector $\vec{x}$ of size $\frac{|\vec{u}|}{2}$ where we sum each two consecutive elements of $\vec{u}$. Now we can see that the $(k-1)$-th element of $\vec{r}$ can be calculated by doing the same process as (i), this time, on $\vec{x}$ (see rows 2 and 3 in RSR++ visualization in Figure 4). We continue this process by repeating steps (i) and (ii) to build all $k$ elements of $\vec{r}$.

In this section, we focus on verifying the theoretical foundation of our proposed methods, RSR and RSR++, using native implementations in C++. We selected C++ as the implementation language because, being a compiled language, it is less susceptible to runtime noise, allowing inference time to better reflect the exact time complexity of the methods, relatively. As baselines, we also implemented Standard Matrix Multiplication, which we will refer to as Standard throughout this section.

For input, we assume access to a binary weight matrix $B^{n\times n}$ where $2^{11}\leq n\leq 2^{16}$ during preprocessing. During inference time, given input vector $\vec{v}$ of size $n$, our aim is to compute $v\cdot B$. with both $B$ and $\vec{v}$ values generated randomly. These parameters align with typical size of weight matrices and input vector patterns commonly encountered in deep learning and large language model (LLM) inference tasks. We then measure inference time across varying values of $n$ for each method. We also utilized the optimal value of $k$ for each $n$, which will be discussed in detail in the next section.

Figure 5 presents the inference times for these methods. Notably, our proposed algorithms achieve up to a $29$x speedup over the Standard baseline when $n=2^{16}$, representing a substantial improvement that could significantly impact DNN and LLM architectures, enabling them to utilize binary and ternary weights more effectively.

Figure 6 illustrates the difference between only RSR++ and RSR for in the inference time. In this figure, we can see up to $25$% improvement when using RSR++ in the Step 2 of inference time (See Algorithm 3). The improvement percentage is computed using the formula $\frac{T({\tt RSR})-T({\tt RSR++})}{T({\tt RSR})}\times 100$.

As outlined in Section 5.2.2, the optimal value of $k$ for RSR (reps. RSR++) is determined through a binary search over the interval $[0,\log(\frac{n}{\log(n)})]$ (resp. $[0,\log(n)]$). Once the optimal $k$ is identified for each $n$, it is then applied in the subsequent experiments. Figures 7(a) 7(b), show the running time with respect to different $k$ values. By increasing $n$, the optimum value of $k$ also increases.

In this section, we investigate the practical performance improvements achieved by RSR, extending beyond native simulation environments. To this end, we use NumPy, a state-of-the-art library for matrix multiplication upon which numerous high-performance computing libraries are based. We implement RSR in Python using NumPy’s built-in functionalities and compare the inference time and memory consumption between RSR and NumPy’s np.dot() function (referred to as NumPy).

• •

  Our RSR++ algorithm shows up to a $29$x improvement in inference time and a $6$x reduction in memory usage, as observed in our experiments. Note that following our complexity analyses, these numbers get even larger by the size of the network layers. This significant improvement indicates that the algorithm can support larger hidden layers and matrices to enhance model accuracy while requiring fewer computational resources.

• •

  The efficiency of our algorithm leads to reduced energy consumption, further contributing to the model’s practicality and sustainability.

• •

  Our approach enables the deployment of more advanced models on ordinary devices, such as personal computers, with limited memory and processing power.

• •

  By processing one column block at a time, the algorithm only requires memory equivalent to the size of a single block, which is significantly less than that needed for full matrix multiplication.

• •

  Last but not least, our method is deployable on top of any of the current or future ternary models, without requiring any fine-tuning or re-training of the LLM or DNN. Once a model is trained, our preprocessing step can be applied, allowing the preprocessed model to be utilized at any future point.

• •

  Our algorithms can be parallelized by executing block productions concurrently, as each block production in step 2 of the inference time is independent. Additionally, further parallelization can be achieved by applying block production simultaneously on different columns of the binary-encoded matrix $Bin_{[k]}$. However, RSR++ is not easily parallelizable due to the sequential nature of operations applied on the vector $\vec{u}$ in step 2. Future work could explore leveraging specific hardware instructions, such as blocking, permutations, segmentations, and block productions, to further optimize these steps at the hardware level.

• •

  The preprocessing time of our algorithms is $O(n^{2})$, and this step must be applied to all weight matrices of an LLM or DNN. Although this process can be time-intensive, it is performed only once for a trained model.

Our proposed algorithms are applicable to any binary or ternary matrices of size $n\times m$. Additionally, they can be extended to any $q$-bit quantized matrix or quantized model. The approach begins by expressing the $q$-bit matrix as the sum of $2^{q-2}$ binary matrices. Each of these binary matrices is processed individually, and the results are combined at the end. The procedure for decomposing the $q$-bit matrix into $2^{q-2}$ binary matrices can be a result of applying Lemma 3.1 recursively on the matrix.