Learning Multiplication-free Linear Transformations

TL;DR

This paper introduces novel dictionary learning algorithms that produce structured, multiplication-free or multiplication-efficient linear transformations, significantly reducing computational complexity for sparse representations.

Contribution

It presents new algorithms for learning structured dictionaries based on factorizations into simple, numerically efficient transformations with closed-form solutions.

Findings

Effective in image data compression

Outperforms traditional transforms in computational efficiency

Achieves comparable or better sparse representation quality

Abstract

In this paper, we propose several dictionary learning algorithms for sparse representations that also impose specific structures on the learned dictionaries such that they are numerically efficient to use: reduced number of addition/multiplications and even avoiding multiplications altogether. We base our work on factorizations of the dictionary in highly structured basic building blocks (binary orthonormal, scaling and shear transformations) for which we can write closed-form solutions to the optimization problems that we consider. We show the effectiveness of our methods on image data where we can compare against well-known numerically efficient transforms such as the fast Fourier and the fast discrete cosine transforms.

Equations

$$\underset{\mathbf{D},\mathbf{X}}{\mathrm{minimize}}$$

$$\text{diag}({\mathbf{D}}^T\mathbf{D})=1_{n\times 1},$$

$$\Vert \text{vec}(\mathbf{X}){\Vert }_0\le sN,$$

$${\mathbf{R}}_{ij}=\left[\begin{array}{ccccc}{\mathbf{I}}_{i-1}& & & & \\ & a& & c& \\ & & {\mathbf{I}}_{j-i-1}& & \\ & b& & d& \\ & & & & {\mathbf{I}}_{n-j}\end{array}\right],\tilde{\mathbf{R}}=\left[\begin{array}{cc}a& c\\ b& d\end{array}\right].$$

$$\mathbf{R}=\prod _{k=1}^m{\mathbf{R}}_{i_k{j}_k}={\mathbf{R}}_{i_m{j}_m}\dots {\mathbf{R}}_{i_1{j}_1}.$$

$$\mathcal{R}_{p}=\left\{x\ \Big{|}\ x=\sum_{t=1}^{p}s_{t}2^{v_{t}};\ s_{t}\in\{\pm 1\},\ v_{t}\in\mathbb{Z}\right\},$$

$${\mathcal{G}}_1=\frac{1}{\sqrt{2}}\{\left[\begin{array}{cc}-1& 1\\ 1& 1\end{array}\right],\left[\begin{array}{cc}1& 1\\ -1& 1\end{array}\right],\left[\begin{array}{cc}1& -1\\ 1& 1\end{array}\right],\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right],$$

$$\left[\begin{array}{cc}1& -1\\ -1& -1\end{array}\right],\left[\begin{array}{cc}-1& -1\\ 1& -1\end{array}\right],\left[\begin{array}{cc}-1& 1\\ -1& -1\end{array}\right],\left[\begin{array}{cc}-1& -1\\ -1& 1\end{array}\right]\},$$

$$\text{and }{\mathcal{G}}_2=\{\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right],\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right],\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],\left[\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right],$$

$$\left[\begin{array}{cc}0& -1\\ -1& 0\end{array}\right],\left[\begin{array}{cc}-1& 0\\ 0& -1\end{array}\right],\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right],\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\}.$$

$$\Vert \mathbf{Y}-{\mathbf{B}}_{ij}\mathbf{X}{\Vert }_F^2=\Vert \mathbf{Y}{\Vert }_F^2+\Vert \mathbf{X}{\Vert }_F^2-2\text{tr}(\mathbf{Z})+C_{ij}^{(t)},$$

$$C_{ij}^{(1)}=$$

$$C_{ij}^{(3)}=$$

$$C_{ij}^{(5)}=$$

$$C_{ij}^{(7)}=$$

$$C_{ij}^{(9)}=$$

$$C_{ij}^{(11)}=$$

$$C_{ij}^{(14)}=$$

$$\mathbf{Z}={\mathbf{YX}}^T\text{ with entries }{Z}_{ij}={\mathbf{y}}_i^T{\mathbf{x}}_j,$$

$$(i^{\star },j^{\star },t^{\star })=\underset{t,i<j}{\mathrm{argmin}}{C}_{ij}^{(t)},$$

$$\mathbf{B}=\prod _{k=1}^m{\mathbf{B}}_{i_k{j}_k}={\mathbf{B}}_{i_m{j}_m}\dots {\mathbf{B}}_{i_1{j}_1}.$$

$$\Vert \mathbf{Y}-\mathbf{BX}{\Vert }_F^2=\Vert {\mathbf{Y}}_k-{\mathbf{B}}_{i_k{j}_k}{\mathbf{X}}_k{\Vert }_F^2,$$

$$\Vert \mathbf{Y}-{\mathbf{O}}_{ij}\mathbf{X}{\Vert }_F^2=\Vert \mathbf{Y}{\Vert }_F^2+\Vert \mathbf{X}{\Vert }_F^2-2\text{tr}(\mathbf{Z})+H_{ij}^{(t)},$$

$$\mathbf{W}=\mathbf{X}{\mathbf{X}}^T\text{ with entries }{W}_{ij}={\mathbf{x}}_i^T{\mathbf{x}}_j.$$

$$\mathbf{M}=\prod _{\ell =1}^q{\bar{\mathbf{M}}}_{\ell },{\bar{\mathbf{M}}}_{\ell }=\prod _{k=1}^{n/2}{\mathbf{B}}_{i_k{j}_k}^{(\ell )},$$

$$\Vert \mathbf{Y}-{\bar{\mathbf{M}}}_1\mathbf{X}{\Vert }_F^2=\Vert \mathbf{Y}{\Vert }_F^2+\Vert \mathbf{X}{\Vert }_F^2-2\text{tr}(\mathbf{Z})+\sum _{k=1}^{n/2}{C}_{i_k{j}_k}^{(t_k)},$$

$${\mathcal{G}}_3=\left\{\left[\begin{array}{cc}1& 0\\ b& 1\end{array}\right],\left[\begin{array}{cc}1& c\\ 0& 1\end{array}\right]\right\},b,c\in \mathbb{R},$$

$$\min \Vert \mathbf{Y}-{\mathbf{S}}_{ij}\mathbf{X}{\Vert }_F^2=\Vert \mathbf{Y}{\Vert }_F^2+\Vert \mathbf{X}{\Vert }_F^2-2\text{tr}(\mathbf{Z})+D_{ij}^{(t)},$$

$$D_{ij}^{(1)}=-(Z_{ji}-W_{ij}{)}^2{W}_{ii}^{-1},D_{ij}^{(2)}=-(Z_{ij}-W_{ji}{)}^2{W}_{jj}^{-1},$$

$$b^{\star }=(Z_{ji}-W_{ij})W_{ii}^{-1}\text{ and }{c}^{\star }=(Z_{ij}-W_{ji})W_{jj}^{-1},$$

$$\Vert$$

$$=$$

$$b_k^{\star }=\frac{{\mathbf{f}}_k^T({\mathbf{x}}_{i_k}\otimes {\mathbf{a}}_{j_k})}{\Vert {\mathbf{x}}_{i_k}{\Vert }_2^2\Vert {\mathbf{a}}_{j_k}{\Vert }_2^2}\text{ and }{c}_k^{\star }=\frac{{\mathbf{f}}_k^T({\mathbf{x}}_{j_k}\otimes {\mathbf{a}}_{i_k})}{\Vert {\mathbf{x}}_{j_k}{\Vert }_2^2\Vert {\mathbf{a}}_{i_k}{\Vert }_2^2},$$

$$\min \Vert \mathbf{Y}-{\mathbf{A}}_k{\mathbf{S}}_{i_k{j}_k}{\mathbf{X}}_k{\Vert }_F^2=\Vert {\mathbf{f}}_k{\Vert }_2^2-E_{ij}^{(t)},\text{ with}$$

$$E_{ij}^{(1)}=\frac{({\mathbf{f}}_k^T({\mathbf{x}}_{i_k}\otimes {\mathbf{a}}_{j_k}){)}^2}{\Vert {\mathbf{x}}_{i_k}{\Vert }_2^2\Vert {\mathbf{a}}_{j_k}{\Vert }_2^2},E_{ij}^{(2)}=\frac{({\mathbf{f}}_k^T({\mathbf{x}}_{j_k}\otimes {\mathbf{a}}_{i_k}){)}^2}{\Vert {\mathbf{x}}_{j_k}{\Vert }_2^2\Vert {\mathbf{a}}_{i_k}{\Vert }_2^2},$$

Full text

I Introduction

In many situations, the success of theoretical concepts in signal processing applications depends on there existing an accompanying algorithmic implementation that is numerically efficient, e.g., Fourier analysis and the fast Fourier transform (FFT) or wavelet theory and the fast wavelet transform (FWT). Unfortunately, in a machine learning scenario where linear transformations are learned they do not exhibit in general advantageous numerical properties, as do the examples just mentioned, unless we explicitly search for such solutions.

In this paper, we propose solutions to the dictionary learning problem [3] which construct linear transformations that have a series of desirable numerical properties while still sparsely representing the training data we supply. Our goal is to build these dictionaries $\mathbf{D}$ such that matrix-vector multiplications $\mathbf{Dx}$ have complexity $O(n\log n)$ or $O(n)$ while we also focus on investigating ways in which the number of multiplication operations can be reduced or completely avoided.

There has been significant work in the literature to learn structured dictionaries that have controllable numerical complexity. One of the earlier attempts is to build a double sparse model [4] where the components of the dictionary are sparse linear combinations from a well-known transform that has a numerically efficient implementation. A recent paper [2] shows how to extend this model and also learn the numerically efficient transformation together with the sparse linear combinations. Other works focus on constructing dictionaries based on Kronecker products [5], circulant (and union of circulants) [6, 7] or convolutional [8] structures, or square transformations that are factored by few Householder reflectors [9], Givens rotations [10] and their generalization [1].

To our knowledge, the dictionary learning community has not investigated the possibility of constructing multiplication free linear transformations. This task has been well studied by the signal and image processing communities where integer-to-integer transformations (also called integer mappings) perform only addition and bit shift operations and therefore are essential for lossless compression. Fast multiplierless approximations of the discrete cosine transform based on a prototype method [11], a lattice structure [12], the lifting scheme [13], an integer [14] and an approximate multiplier-less [15] fast Fourier transform were developed first. Then [16] introduced a general framework to build integer mappings from any linear transformation based on factorizations of (triangular and row) elementary reversible matrices and then showcases the framework on the discrete Fourier, cosine and wavelet transformations. Another general framework based on the general S transform is given in [17]. One image processing application is for the design of an integer color transform [18].

In this paper, we will combine the benefits from both worlds: in the style of dictionary learning, we will learn a numerically efficient transformation from a training dataset that directly has an imposed structure to reduce or eliminate multiplication operations, in the fashion of integer mappings.

The paper is structured as follows: in Section II we discuss ways to measure computational complexity, in Section III we briefly describe the dictionary learning problem and our computational design goals, then in Sections IV and V we develop the proposed learning procedures and finally in Section VI we show experimental results with image data where we compare to the discrete cosine transform.

II A note on computational complexity

Given the number of computational platforms available today and their sophistication, our purpose is not to provide an exhaustive, detailed discussion of the subject but to give arguments that multiplication-free algorithms are relevant.

In most scenarios the computational complexity accounts for all the operations performed by the system, i.e., we count together mathematical operations like additions, multiplications etc. When considering modern computing systems this choice is a natural one: these mathematical operations take approximately the same time sophisticated hardware. For example, numerical simulation performed with an Intel i$7^{\text{\tiny{\textcopyright}}}$ processor shows that integer multiplication is on average approximately only 10% slower than integer addition (running Linux, using the gcc with the –O3 flag, the program performs the operations on random integer operands in two arrays and the results are stored in a third). Historically, this was not the case. In the past, computer scientists have made several efforts to reduce the number of multiplications in their algorithmic implementations in favor of performing more addition operations. A classic example is the multiplication of two complex numbers which can be done in two ways: the first takes four multiplications and two additions while the second one has three multiplications and five additions, i.e., it was computationally convenient to replace one multiplication by three addition operations. Even so, modern computing systems still perform integer addition faster than integer multiplication in general (one clock cycle versus three to ten clock cycles depending on the particular processor)111Intel 64 and IA-32 Architectures Optimization Reference Manual.

From an algorithmic perspective, for numbers represented using $n$ bits, it is well understood that integer addition has complexity $O(n)$ while the best asymptotic bound $O(n\log n\log\log n)$ for integer multiplication, given for realistically reasonably large $n$, is achieved by the Schonhage-Strassen algorithm [19]. Other, asymptotically less efficient approaches, include Karatsuba’s algorithm [20] and, its generalization, the Toom-Cook algorithm [21, Section 9.5] – both use techniques similar to the previously described trick of replacing an intermediate multiplication operation with several additions achieve complexity $n^{\log 3}$ and $O\left(n^{\frac{\log 5}{\log 3}}\right)$, respectively. In terms of hardware, the modern multiplier architectures use the Baugh-Wooley algorithm [22] or Wallace trees [23], for example. These methods reduce the performance gap with binary addition (as previously observed experimentally) at the cost of increasing the complexity of the circuitry.

Aside from the execution time, there are several other important complexity measures, like power consumption and circuitry size, especially when considering some custom or embedded computational platforms where low size, weight, power and cost (SWaP-C) solutions are preferred. An $n$ bit full-adder needs $5n$ logic gates: one OR, two AND and two XOR gates per bit. In the case of binary multiplication, for example, the relatively simple sequential $n$ bit array multiplier needs $31n$ gates: the $n$ and $2n$ bit registers consist of $15n$ gates, the ALU contains an adder and two multiplexers which consist of $16n$ gates ($5n$ gates for the adder, $4n$ gates for the $2\times 1$ mux and $7n$ gates for the $4\times 1$ mux).

Application-specific integrated circuits (ASICS) are circuits that are designed to perform only one (or a small set) of tasks, unlike CPUs. Field-programmable gate arrays (FPGAs) are a computational platform that belongs to the ASICS class. Addition operations in FPGAs are generally performed using look-up tables while for the multiplications some specialized extra components are needed (like DSP slices). Here, a frequently used performance indicator is the power-delay product (the product between the energy consumption and the input-output delay of a circuit). For a very popular FPGA computer-aided design tools for arithmetic code generation, the $\text{Xilinx}^{\text{\tiny{\textcopyright}}}$ IP Core Generator, with 32-bit operands the addition operation has a power-delay product of 0.67 nJ (see Table 4.1 of [24]) while the 15-bit multiplication operation has a power-delay product of 2.21 nJ (see Table 4.3 of [24]).

The same power-delay product values for microprocessors at 45nm are given in [25] to be: 0.1pJ and 3pJ for 32-bit integer addition and multiplication, respectively, and 0.9pJ and 4pJ for 32-bit floating point addition and multiplication, respectively. These numbers do not take into account memory access latency and power consumption (which according to [25] are also a major contributor to the overall power consumption).

III Preliminaries

In this section, we briefly describe the general dictionary learning problem and then proceed to list some computationally-desirable properties of the learned dictionary. We also describe some basic matrix building blocks that we will use in this paper to reach the desirable properties listed.

III-A The dictionary learning problem

Given a $N$ point dataset $\mathbf{Y}\in\mathbb{R}^{n\times N}$ and the average sparsity $s\in\{1,\dots,n-1\}$, the dictionary learning problem can be stated in the optimization language as

[TABLE]

where $\mathbf{D}\in\mathbb{R}^{n\times n}$, which always has unit $\ell_{2}$ columns, is called the dictionary (in general $\mathbf{D}$ can be overcomplete, but in this paper we consider square dictionaries which we also call transforms [26]) and $\mathbf{X}\in\mathbb{R}^{n\times N}$ has the sparse representations of all data points in the dictionary $\mathbf{D}$. The $\ell_{0}$ pseudo-norm constraint, which counts the number of non-zero in the matrix $\mathbf{X}$, ensures that on average each data point is represented using $s$ columns (also called atoms) from $\mathbf{D}$. The dictionary learning problem is hard in general, so most optimization techniques lead to local minima points of (1) by a process of alternating minimization: keep $\mathbf{D}$ fixed and create $\mathbf{X}$ and then vice-versa. In this paper, we deploy the same alternating technique and we focus on constructing the dictionary $\mathbf{D}$ such that it has some specific computational properties. We do not focus on how to construct the sparse representations $\mathbf{X}$, i.e., we will use the appropriate, well-established algorithms from the literature to build $\mathbf{X}$ [27, Chapter 1].

For a recent, detailed description of the dictionary learning problem and some of its solutions the reader is encouraged to check [27, Chapters 2 and 3].

III-B The basic building block

Based on the work in [1], we revise the $n\times n$ R-transform:

[TABLE]

These matrices can be viewed as perturbations of the identity matrix: $\mathbf{R}_{ij}$ has zero entries everywhere except for its diagonal (with entries $a$ and $d$ in positions $(i,i)$ and $(j,j)$, respectively, and the rest with value one) and the only two off-diagonal entries $c$ and $b$, on positions $(i,j)$ and $(j,i)$ respectively. The subscripts of the R-transform define the rows on which the non-trivial values are stored. As it is convenient in many situations to reference the unique part of $\mathbf{R}_{ij}$ separately we denote it by $\mathbf{\tilde{R}}$. Through this paper we will use the same template for the transformations we propose: $n\times n$ identity except for 2 coordinates where we will perform a carefully chosen $2\times 2$ calculation. We define transformations that are products of $m$ basic building blocks, like

[TABLE]

Analogously to (2), we denote the unique $2\times 2$ part of each $\mathbf{R}_{i_{k}j_{k}}$ by $\mathbf{\tilde{R}}_{k}$.

Given $\mathbf{x}\in\mathbb{R}^{n}$ the matrix-vector multiplication $\mathbf{R}_{i_{k}j_{k}}\mathbf{x}$ takes four multiplications and two addtions. To avoid these multiplications, similarly to fixed point number representations, we use the sums of powers of two (SOPOT) set:

[TABLE]

where the parameter $p$ establishes the precision of the entries. By convention when $p=\infty$ we use the precision of the working data type (double floating point in our case). For given $p$, we provide in Algorithm 1 an iterative greedy procedure to compute the representation of any real scalar input $x$ in the set $\mathcal{R}_{p}$. We use this set to represent our transformations $\mathbf{\tilde{R}}$ which we now denote $\mathbf{\tilde{R}}_{k,p}\in\mathcal{R}_{p}^{2\times 2}$, i.e., we approximately represent each $a,b,c,d$ in the set $\mathcal{R}_{p}$, and we call the overall $\mathbf{R}_{i_{k}j_{k},p}$ an Rp-transform.

Structures like these Rp-transforms are interesting numerically because matrix-vector multiplication takes $4p$ bit shifts and $4p-2$ additions ($4p-4$ to form the four products and 2 to add the results for line $i$ and $j$ respectively). This is because scalar multiplication with $a\in\mathcal{R}_{p}$ takes $p$ bit shifts and $p-1$ additions. In terms of the coding complexity, assuming $8$ bits are used to store each $v_{t}$ in (4), storing $\mathbf{R}_{i_{k}j_{k},p}$ takes approximately $36p+2\log_{2}n$ bits ($2\log_{2}n$ bits to store the indices $i,j$ and $9p$ bits to store each entry in $\mathcal{R}_{p}$).

In this paper, we propose to learn dictionaries which are products of basic transformations like (2), while we also impose some additional constraints, e.g., orthogonality or some specific arithmetic structure for the non-zero entries.

III-C The computational properties for the dictionary

In this section we define and present properties of the basic building blocks we consider for numerically efficient factorizations that will allow us to achieve our design goals.

Our goal is to construct a dictionary $\mathbf{D}\in\mathbb{R}^{n\times n}$ that display properties such as:

1. P1.

The computational complexity (the number of additions, bit shifts and multiplications) of $\mathbf{Dx}$ and $\mathbf{D}^{-1}\mathbf{x}$ for any given $\mathbf{x}\in\mathbb{R}^{n}$ is controllable, preferably $O(n\log n)$; 2. P2.

If $\mathbf{x}\in\mathcal{R}_{p}^{n}$ then $\mathbf{Dx}\in\mathcal{R}_{p\prime}^{n}$, i.e., if $\mathbf{x}$ has a fixed point representation then so does $\mathbf{Dx}$; if $\mathbf{Dx}$ is calculated exactly in $\mathcal{R}_{p}^{n}$ then $\mathbf{D}^{-1}\mathbf{x}$ can also be calculated exactly in $\mathcal{R}_{p\prime}^{n}$ with $p\prime\neq p$, i.e., if $\mathbf{D}$ has a fixed representation then so does the inverse $\mathbf{D}^{-1}$; 3. P3.

$\mathbf{D}$ is exactly reversible, i.e., $\mathbf{D}^{-1}\mathbf{D}=\mathbf{I}$, when $\mathbf{D}$ has a fixed point representation; 4. P4.

Reduce, and ideally, eliminate multiplication operations for $\mathbf{Dx}$ and $\mathbf{D}^{-1}\mathbf{x}$;

In the following sections, we will distinguish between orthonormal (Section IV) and general (Section V) dictionary learning procedures and discuss how the transformations we learn achieve some or all of these desirable properties.

IV The orthonormal case

In this section, we propose two orthonormal dictionary learning algorithms: one with a reduced number of multiplications and one that avoids completely such operations.

IV-A Numerically efficient orthogonal transforms: Bm–DLA

We define the two sets of orthonormal binary $2\times 2$ matrices:

[TABLE]

[TABLE]

We define an orthogonal B-transform denoted $\mathbf{B}_{ij}$ as a constraint R-transform in (2) where we have the non-trivial, non-zero $2\times 2$ part positioned at indices $i$ and $j$, which we denote $\mathbf{\tilde{B}}_{ij}$, defined as one of the sixteen options, i.e., $\mathbf{\tilde{B}}_{ij}\in\{\mathcal{G}_{1}\cup\mathcal{G}_{2}\}$. Notice that given $\mathbf{A},\mathbf{B}\in\mathcal{G}_{2}$ we have that $\mathbf{AB}\in\mathcal{G}_{2}$, $\mathbf{BA}\in\mathcal{G}_{2}$ (and in fact $\mathcal{G}_{1}\cup\mathcal{G}_{2}$ has a group structure) and given $\mathbf{C}\in\mathcal{G}_{1}$ we have that $\mathbf{AC}\in\mathcal{G}_{1}$ and $\mathbf{CA}\in\mathcal{G}_{1}$. The factor $2^{-\frac{1}{2}}$ is there to keep the B-transforms orthonormal (orthogonal and with columns normalized in $\ell_{2}$).

Structures like B-transforms are useful because matrix-vector multiplication takes four operations: two additions and two multiplications (both by $2^{-\frac{1}{2}}$). The coding complexity of storing $\mathbf{B}_{ij}\in\mathbb{R}^{n\times n}$ is approximately $4+2\log_{2}n$ bits (the first term encodes the choice in $\mathcal{G}_{1}\cup\mathcal{G}_{2}$ while the second encodes the two indices $i$ and $j$). We could avoid the multiplications by approximating $2^{-\frac{1}{2}}$ in $\mathcal{R}_{p}$ but we would lose orthogonality.

We are interested in solving optimization problems that consider one B-transform as a dictionary:

[TABLE]

where we have used the definition of the Frobenius norm $\|\mathbf{A}\|_{F}^{2}=\text{tr}(\mathbf{A}^{T}\mathbf{A})$ and the index $t=1,\dots,15,$ runs through the possible variants $\mathbf{\tilde{B}}_{ij}$ in (5), (6). Therefore, we have for the order of the transformations in (5), (6):

[TABLE]

with $Z_{sd}=Z_{ii}+Z_{jj},Z_{so}=Z_{ij}+Z_{ji},Z_{do}=Z_{ij}-Z_{ji}$, the constants $c_{1}=2+\sqrt{2},\ c_{2}=2-\sqrt{2}$ and we define:

[TABLE]

where $\mathbf{y}_{i}^{T}$ and $\mathbf{x}_{i}^{T}$ are the $i^{\text{th}}$ rows of $\mathbf{Y}$ and $\mathbf{X}$, respectively.

The $\mathbf{B}_{ij}$ that minimizes (7) is given by

[TABLE]

for the $C_{ij}^{(t)}$ in (8) with $t=1,\dots,15,$ and $j=1,\dots,n-1$. The total computational complexity to find the minimizer of (10) is: $2n^{2}N$ operations to construct $\mathbf{Z}$, which dominates the computational complexity; $O(n^{2})$ operations to compute all the fifteen $C_{ij}^{(t)}$ for all $\frac{n(n-1)}{2}$ distinct pairs $(i,j)$ with $i<j$ and solve (10); and $2n$ operations to compute $2\text{tr}(\mathbf{Z})$.

We propose a method to learn orthogonal dictionaries that are factorized as a product of $m$ B-transforms. Therefore we propose the following structure for our learned dictionary:

[TABLE]

With this choice, the dictionary learning objective function for a single transformation indexed $k$ is:

[TABLE]

where $\mathbf{Y}_{k}=\prod_{q=k+1}^{m}\mathbf{B}_{i_{q}j_{q}}^{T}\mathbf{Y}$ and $\mathbf{X}_{k}=\prod_{q=1}^{k-1}\mathbf{B}_{i_{q}j_{q}}\mathbf{X}$. In this development we have used the fact that orthonormal transformations are invariant in the Frobenius norm, i.e., $\|\mathbf{QY}\|_{F}=\|\mathbf{Q}^{T}\mathbf{Y}\|_{F}=\|\mathbf{Y}\|_{F}$ for any orthonormal $\mathbf{Q}\in\mathbb{R}^{n\times n}$. Notice that we have reduced the objective function to the form in (7). Therefore, we propose an efficient iterative process that updates a single $\mathbf{B}_{i_{k}j_{k}}$ at a time while keeping the others fixed until all components are optimized.

We describe the full learning procedure in Algorithm 2. This algorithm updates iteratively each B-transform in the composition of the dictionary $\mathbf{B}$ and then the sparse representations $\mathbf{X}$. Since each step is solved exactly to optimality, the algorithm monotonically converges overall to a local minimum point or stops early in the maximum $K$ iterative steps. In the description of the algorithms, we have used $\mathcal{T}_{s}()$, an operator that given a matrix sets to zero, in each column separately, all entries except the highest $s$ in magnitude.

Remark 1 (Working with limited precision). In order to design exactly invertible orthogonal linear transformations when using data with fixed bit representation, lifting schemes [28, 29] were introduced in the past: $\begin{bmatrix}c&-s\\ s&c\end{bmatrix}=\begin{bmatrix}1&\frac{c-1}{s}\\ 0&1\end{bmatrix}\begin{bmatrix}1&0\\ s&1\end{bmatrix}\begin{bmatrix}1&\frac{c-1}{s}\\ 0&1\end{bmatrix},\begin{bmatrix}c&s\\ s&-c\end{bmatrix}=\begin{bmatrix}1&\frac{c-1}{s}\\ 0&1\end{bmatrix}\begin{bmatrix}1&0\\ s&1\end{bmatrix}\begin{bmatrix}1&-\frac{c-1}{s}\\ 0&-1\end{bmatrix}$. The proposed B-transforms can naturally be implemented with these schemes, as they are particular $2\times 2$ orthonormal matrices. For the transformations in $\mathcal{G}_{1}$ we have that $\frac{c-1}{s}\in\{\pm(1-\sqrt{2}),\pm(1+\sqrt{2})\}$ and $s\in\left\{\pm 2^{-\frac{1}{2}}\right\}$ while the transformations in $\mathcal{G}_{2}$ are multiplier-less and therefore do not need the lifting scheme representations. Still, note that matrix-vector multiplications with matrices from $\mathcal{G}_{1}$ take two multiplications and two additions, while using lifting schemes representations needs three multiplications and three additions. The elements $\frac{c-1}{s}$ and $s$ in the lifting scheme can be represented in $\mathcal{R}_{p}$ to avoid multiplications altogether, but at the cost of loosing orthogonality.$\hfill\blacksquare$

Remark 2 (Avoiding the normalization factor). The normalization by $2^{-\frac{1}{2}}$ seems to cause complications from a numerically efficiency stand point. Therefore, we could define an O-transform $\mathbf{O}_{ij}\in\mathbb{R}^{n\times n}$ which is achieved for (2) when $a,b,c,d\in\{\pm 1\}$ such that $\mathbf{\tilde{R}}_{ij}$ is orthogonal. These eight transformations use $\sqrt{2}\mathbf{\tilde{B}}_{ij}$ for the structure (5). Solving such least squares problems leads to

[TABLE]

where, for $t=1,\dots,8$, we have defined $H_{ij}^{(t)}=W_{ii}+W_{jj}-2(a-1)Z_{ii}-2(d-1)Z_{jj}-2bZ_{ij}-2cZ_{ji}+2(ab+cd)W_{ij}$. We have used the notation:

[TABLE]

where $\mathbf{x}_{i}^{T}$ is the $i^{\text{th}}$ row of $\mathbf{X}$.

Similarly to (11), based on these simple transformation we can define $\mathbf{O}=\prod_{k=1}^{m}\mathbf{O}_{i_{k}j_{k}}=\mathbf{O}_{i_{m}j_{m}}\dots\mathbf{O}_{i_{1}j_{1}}.$ Because the normalization is entirely avoided, $\mathbf{O}_{ij}$ is no longer orthogonal and therefore $\mathbf{O}$ is not orthogonal. A drawback of this is the fact that the transformations $\mathbf{O}_{i_{k}j_{k}}$ cannot be rearranged as in (12), making the update of an individual transformation while keeping all others fixed more difficult.

But this structure has the advantage of completely avoiding multiplication operations, i.e., matrix-vector multiplications $\mathbf{O}_{ij}\mathbf{x}$ take 2 addition operations. Also, notice that $\mathbf{O}$ has integer entries and that $\det(\mathbf{O})=\pm 2^{m}$ since $\det(\mathbf{O}_{ij})=\pm 2$.

A sufficient condition for the local optimality of $\mathbf{O}_{ij}$ is that $\|\mathbf{x}_{j}\|_{2}-\|\mathbf{y}_{i}-\mathbf{x}_{i}\|_{2}\geq 0,\ \forall\ i\neq j$, i.e., the energy of any row error is less than the energy of all other rows of $\mathbf{X}$.

Proof. Given any rows of $\mathbf{Y}$ and $\mathbf{X}$ there is no O-transform that improves the objective function if $\|\mathbf{y}_{i}-\mathbf{x}_{i}\|_{2}^{2}\leq\|\mathbf{y}_{i}-(\mathbf{x}_{i}\pm\mathbf{x}_{j})\|_{2}^{2},\ \forall\ i\neq j$. Developing this leads to $\|\mathbf{x}_{j}\|_{2}^{2}-|C|\geq 0,\ C=2\mathbf{x}_{i}^{T}\mathbf{x}_{j}-2\mathbf{y}_{i}^{T}\mathbf{x}_{j}$. By the Cauchy-Schwartz inequality we have that $-|C|\geq-\|\mathbf{x}_{j}\|_{2}\|\mathbf{y}_{i}-\mathbf{x}_{i}\|_{2}$ and therefore $\|\mathbf{x}_{j}\|_{2}^{2}-|C|\geq\|\mathbf{x}_{j}\|_{2}^{2}-\|\mathbf{x}_{j}\|_{2}\|\mathbf{y}_{i}-\mathbf{x}_{i}\|_{2}\geq 0$. $\hfill\blacksquare$

The matrix-vector multiplication operation with the whole $\mathbf{B}$ transformation takes $2m$ additions and $2m$ multiplications with the same constant value $2^{-\frac{1}{2}}$.

Related to our previously stated desired computational properties we have that: P1 is achieved by taking the number of B-transforms in $\mathbf{B}$ to be $m\sim O(n\log n)$ and since $\mathbf{B}^{-1}=\mathbf{B}^{T}$ we have that the inverse transformation enjoys the same computational complexity (the inverse of any B-transforms is itself a B-transform); P2 can be achieved by using the lifting schemes and representing $\frac{c-1}{s}$ and $s$ in $\mathcal{R}_{p}$ but notice that with the fixed point representations we no longer have an orthogonal transformation, we denote $\mathbf{B}_{p}$ the transformation $\mathbf{B}$ in the lifting scheme with elements in $\mathcal{R}_{p}$; using the representation explained for P2 we also cover the requirement P3; regarding P4, notice that multiplication with each B-transform takes two multiplications as compared to three in the general lifting scheme.

IV-B Multiplication-free orthogonal transforms: Mq–DLA

B-transforms are numerically efficient structures, although they do involve multiplication operations. In this section we explore ways to reduce the numerical complexity even further. Consider the following structure:

[TABLE]

with $\bigcup_{k=1}^{n/2}(i_{k},j_{k})=\{1,\dots,n\}\text{ and }\bigcap_{k=1}^{n/2}(i_{k},j_{k})=\emptyset$. This transformation is made up of $q$ stages. At each stage there are $\frac{n}{2}$ B-transforms that are chosen such that the coordinates $(i_{k},j_{k})$ are a partition of $\{1,\dots,n\}$. Consider now the objective function for a single block $\mathbf{\bar{M}}_{1}$

[TABLE]

such that the indices $(i_{k},j_{k})$ obey the constraints in (15). Optimizing the expression in (16), i.e., minimizing the summation term by finding the best parings of the indices, is equivalent to the weighted maximum matching algorithm [30] (of maximum-cardinality matchings) on the graph with $n$ nodes and with edges $C_{i_{k}j_{k}}^{(0)}=\underset{t_{k}=1,\dots,8}{\min}-\ C_{i_{k}j_{k}}^{(t_{k})}$ (the minus sign is set because we want to minimize the quantity). Because all $\mathbf{\bar{M}}_{\ell}$ are orthonormal, the manipulations in (12) hold and therefore each $\mathbf{\bar{M}}_{\ell}$ can be updated while the others are fixed. We have found in our experimental settings that these iterative steps do not significantly improve the solution reached and therefore the algorithm builds the dictionary $\mathbf{M}$ in a single iteration, i.e., we set $K=1$. This also highlights the importance of the initialization for $\mathbf{X}$, which is done again by the singular value decomposition.

The full procedure is shown in Algorithm 3. The main difference with the previously introduced Bm–DLA is that at each step we update $\frac{n}{2}$ B-transforms simultaneously, not just one and none of these transforms use the same coordinates. In this fashion, the constant $2^{-\frac{1}{2}}$ factors out while keeping the $\mathbf{M}$ transformation orthonormal.

The transformation in (15) can be equivalently written as $\mathbf{M}=2^{-\frac{q}{2}}\prod_{\ell=1}^{q}\mathbf{M}_{\ell}$, where all $\mathbf{M}_{\ell}\in\mathbb{R}^{n\times n}$ are orthogonal sparse matrices with elements in $\{0,\pm 1\}$, the diagonal only in $\{\pm 1\}$ and two non-zero entries per each row and column. Therefore, matrix-vector multiplication with one $\mathbf{M}_{\ell}$ takes $n$ additions and one multiplication and as such, matrix-vector multiplication with the whole $\mathbf{M}$-transform takes $nq$ additions and $n$ multiplications (or bit shifts if $q$ is even). From this description of $\mathbf{M}$ it is easy to see that the coding complexity is approximately $\frac{q}{\ln 2}(n\ln n-n+1)$ bits: the cost of encoding $q$ partitions of the indices, i.e., $q\sum_{i=1}^{n}\log_{2}i=q\log_{2}n!\approx\frac{q}{\ln 2}(n\ln n-n+1)$ by Stirling’s approximation. The constant $q$ is encoded implicitly as the number of partitions.

The significant benefit of transforms like (15) is that they avoid multiplications altogether. Furthermore, notice that each stage can be completely parallelized (since operations are done on distinct indices). The disadvantage, especially when compared to Bm–DLA, is that we force the transform to use all available indices and this constraint is sub-optimal in general and inferior to the choice (10) made in Bm–DLA. Therefore, we expect Mq–DLA to perform worse than Bm–DLA in terms of the objective function value for the same number of basic transformations, i.e., $m=\frac{nq}{2}$.

The computational complexity of Mq–DLA is dominated by the computation of $\mathbf{Z}$ which takes $O(n^{2}N)$ operations and by the overall iterative process which takes a total of $O(qn^{3})$ ($q$ times we have to perform the partitioning of the indices by the maximum matching algorithm).

Remark 3 (On the computational complexity of finding the best partition of indices). The weighted maximum matching algorithm has complexity $O(n^{3})$ which might be prohibitively large in some learning situations – especially as compared to Bm–DLA which has a $O(n^{2})$ complexity per iteration. An alternative is to use a sub-optimal, greedy, approach to build the indices partition. Consider a procedure that builds the partition in two iterative steps: compute the $\mathbf{B}^{(\ell)}_{i_{k}j_{k}}$ by $(i_{k},j_{k},t_{k})=\underset{t,i<j;\ i,j\notin\mathcal{S}}{\arg\min}\ -2\text{tr}(\mathbf{Z}_{\ell})+C_{ij}^{(t)}$ and then update the set $\mathcal{S}\leftarrow\mathcal{S}\cup(i_{k},j_{k})$ for $k=1,\dots,\frac{n}{2}$ starting from $\mathcal{S}=\emptyset$. $\hfill\blacksquare$

Remark 4 (Another strategy for avoiding multiplication operations). Notice that the matrices in (5) are permutations with sign flips of the $2\times 2$ Hadamard matrix $\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\ 1&-1\end{bmatrix}$. We can extend these structures by using the $4\times 4$ Hadamard matrix $\frac{1}{2}\left(\begin{bmatrix}1&1\\ 1&-1\end{bmatrix}\otimes\begin{bmatrix}1&1\\ 1&-1\end{bmatrix}\right)$, whose scaling factor is now simple, i.e., a power of two, and avoids multiplication operations. Unfortunately, operating on more than two coordinates increases the numerical complexity of the learning procedure, i.e., for each pair of four indices, instead of 8 options in (5) we now have $768$ options, equivalent to all possible permutations of the four rows and columns and sign changes ($2\times 2^{4}\times 4!$) and there are ${n\choose 4}\approx\frac{n^{4}}{24}$ such pairs of indices, instead of $\frac{n(n-1)}{2}$ as for Bm–DLA. The total overall cost of one training iteration would therefore be dominated by the computation of the approximately $32n^{4}$ quantities $C_{i_{k}j_{k}r_{k}p_{k}}^{(t_{k})}$ from which the minimum has to be found. Operating on even more coordinates simultaneously seems unreasonable from a computational perspective (in the learning phase).

Still, the benefit is that the matrix-vector multiplication with a single such structure takes $12$ addition and $4$ bit shift operations. Therefore, an algorithm that uses these fundamental building blocks produces a transformation that completely avoids the multiplication operations. Finally, note that there are 5378 possibilities when considering all $3\times 3$ orthonormal matrices (with different scaling factors $\sqrt{2}$, $\sqrt{3}$ and $2$) with entries in $\{0,\pm 1\}$, i.e., if we also allow zero entries – these structure include (5) and (6). Therefore, this structure can be used only in scenarios where the learning time is not fundamentally constrained by time or power considerations.

The advantage is that, as with the other algorithms described in this paper, parallelization is trivial. Furthermore, this approach would combine two of the major benefits of Bm–DLA and Mq–DLA: no calculations of partitions are necessary, i.e., we are not forced to repeatedly use coordinates that do not lead to significant reductions in the objective function, and there are no multiplication operations.

We call this approach B${}_{m}^{\otimes}$–DLA, and it follows the same steps as Bm–DLA but for only $K=1$ number of iterations. For brevity we omit the full description of the algorithm.$\hfill\blacksquare$

Related to our previously stated desired computational properties we have that: P1 is achieved by taking $q\sim O(\log n)$ and since $\mathbf{M}^{-1}=\mathbf{M}^{T}$ we have the same computational benefits for the inverse transformation; P2, P3, and P4 are trivially achieved when $q$ is fixed to be even and therefore the algorithm avoids completely any multiplication operations.

V The general case

In this section, we propose an algorithm to learn general dictionaries which have controllable complexity, including a variant which completely avoids multiplication operations. We begin by discussing the properties of scaling and shear transformations and then propose the learning procedure.

V-A Shear transformations

Consider the set of shear $2\times 2$ matrices

[TABLE]

and define a shear transformation $\mathbf{S}_{ij}\in\mathbb{R}^{n\times n}$ which is achieved for (2) when $\mathbf{\tilde{R}}_{ij}\in\mathcal{G}_{3}$, i.e., $a=1$ and $d=1$ fixed while $b$ and $c$ are free parameters or set to zero, alternatively.

The objective function of our learning problem when the dictionary is a single $\mathbf{S}_{ij}$ now leads to

[TABLE]

for $t\in\{1,2\}$, where $D_{ij}^{(1)}=b^{2}W_{ii}+2b(W_{ij}-Z_{ji})$ and $D_{ij}^{(2)}=c^{2}W_{jj}+2c(W_{ji}-Z_{ij})$. The minima, for $i=1,\dots,n$ and $j=i+1,\dots,n$, are

[TABLE]

and are achieved for the optimum choices

[TABLE]

respectively. Starting with all transformations $\mathbf{S}_{i_{k}j_{k}}$ set to the identity matrix, each one of the transformations $\mathbf{S}_{i_{k}j_{k}}$ is initialized sequentially in this fashion for $k=1,\dots,m$.

Remark 5 (Optimality condition). A necessary and sufficient condition for local optimality is that $D_{ij}^{(1)}=D_{ij}^{(2)}=0$ and therefore $\mathbf{x}_{j}^{T}(\mathbf{y}_{i}-\mathbf{x}_{i})=0,\ \forall\ i\neq j$, i.e., in the spirit of the least squares solution applied row-wise, we have that any error row $\mathbf{\epsilon}_{i}=\mathbf{y}_{i}-\mathbf{x}_{i}$ is orthogonal to all rows $\mathbf{x}_{i}^{T}$ of $\mathbf{X}$.$\hfill\blacksquare$

Now, after the initialization process, each $\mathbf{S}_{i_{k}j_{k}}$ is updated again while all other S-transforms are kept fixed. The objective function develops now to

[TABLE]

where $\mathbf{y}=\text{vec}(\mathbf{Y})$, $\mathbf{B}_{k}=\mathbf{X}_{k}^{T}\odot\mathbf{A}_{k}\in\mathbb{R}^{nN\times n}$, $\odot$ is the Khatri-Rao product, $\mathbf{f}_{k}=\mathbf{y}-\mathbf{B}_{k}\mathbf{1}_{n\times 1}$, $\mathbf{A}_{k}=\prod_{q=k+1}^{m}\mathbf{S}_{i_{q}j_{q}}$, $\mathbf{X}_{k}=\prod_{q=1}^{k-1}\mathbf{S}_{i_{q}j_{q}}\mathbf{X}$, $\mathbf{L}_{i_{k}j_{k}}\in\{b_{k}\mathbf{e}_{j_{k}}\mathbf{e}_{i_{k}}^{T},c_{k}\mathbf{e}_{i_{k}}\mathbf{e}_{j_{k}}^{T}\}$ and $\{\mathbf{e}_{i}\}_{i=1}^{n}$ are the standard basis vectors of $\mathbb{R}^{n}$. We have used that $\text{vec}(\mathbf{ABC})=(\mathbf{C}^{T}\otimes\mathbf{A})\text{vec}(\mathbf{B})$. Notice that $(\mathbf{X}_{k}^{T}\otimes\mathbf{A}_{k})\text{vec}(\mathbf{e}_{j_{k}}\mathbf{e}_{i_{k}}^{T})$ selects the $(j_{k}+(i_{k}-1)n)^{\text{th}}$ column, i.e., $\mathbf{x}_{i_{k}}\otimes\mathbf{a}_{j_{k}}$, while $(\mathbf{X}_{k}^{T}\otimes\mathbf{A}_{k})\text{vec}(\mathbf{e}_{i_{k}}\mathbf{e}_{j_{k}}^{T})$ selects the $(i_{k}+(j_{k}-1)n)^{\text{th}}$ column, i.e., $\mathbf{x}_{j_{k}}\otimes\mathbf{a}_{i_{k}}$.

To minimize the quantity in (21) the optimal choices are

[TABLE]

respectively, and the minimum objective function values are

[TABLE]

[TABLE]

for all $i=1,\dots,n$ and $j=i+1,\dots,n$.

As with the previously introduced structures, shear transformations have good numerical properties, i.e., matrix-vector multiplication $\mathbf{S}_{i_{k}j_{k}}\mathbf{x}$ takes one addition and one multiplication operation. If the coefficients $b$ or $c$ are represented in $\mathcal{R}_{p}$ then $\mathbf{S}_{i_{k}j_{k},p}\mathbf{x}$ takes $p$ bit shifts and $p$ additions. Furthermore, inverses $\mathbf{S}_{i_{k}j_{k}}^{-1}$ shears themselves and easy to compute because $\begin{bmatrix}1&0\\ b&1\end{bmatrix}^{-1}=\begin{bmatrix}1&0\\ -b&1\end{bmatrix}$ and $\begin{bmatrix}1&c\\ 0&1\end{bmatrix}^{-1}=\begin{bmatrix}1&-c\\ 0&1\end{bmatrix}$ and they have the same numerical properties as the direct shear transformations. There is no scaling factor for these inverses since $\det(\mathbf{S}_{i_{k}j_{k}})=1$ always.

To encode one shear transformation we need approximately $1+C+2\log_{2}n$ bits (one bit to encode the choice between the two shears in (17), the constant $C$ is the cost of encoding of $b$ or $c$, say $C=64$ for a double float, while the second term encodes the indices $i$ and $j$).

V-B Scaling transformations

Consider the R-transform in (2) constrained to $b=0,c=0$, $d=1$ and we drop the index $j$ which is now unnecessary to obtain a scaling matrix along a single coordinate

[TABLE]

Similarly to (18), with this scaling as the dictionary, the objective function of our learning problem is now

[TABLE]

where we have denoted $F_{i}=-2Z_{ii}(a^{\star}-1)+W_{ii}((a^{\star})^{2}-1)$ and used the scalar least squares solution $a^{\star}=Z_{ii}W_{ii}^{-1}$ that minimizes $\|\mathbf{y}_{i}-a^{\star}\mathbf{x}_{i}\|_{F}^{2}$. Our goal is to find $\alpha$ such that $|\alpha a^{\star}|=2^{\delta},\ \delta\in\mathbb{Z},$ and $\|\mathbf{y}_{i}-\alpha a^{\star}\mathbf{x}_{i}\|_{F}^{2}$ is minimized. It is necessary to verify that the further scaling $\alpha$ is such that $\|\mathbf{y}_{i}-\mathbf{x}_{i}\|_{F}^{2}-\|\mathbf{y}_{i}-\alpha a^{\star}\mathbf{x}_{i}\|_{F}^{2}\geq 0$, which is obeyed when $-\left|\frac{1-a^{\star}}{a^{\star}}\right|\leq\alpha-1\leq\left|\frac{1-a^{\star}}{a^{\star}}\right|$, i.e., our scaling does not increase the objective function as compared to doing nothing along the $i^{\text{th}}$ coordinate. The left-hand side of the previous inequality is maximized when $\alpha=2^{\delta}|a^{\star}|^{-1}$ is closest to one and therefore

[TABLE]

As such, the minimizer of (26) that is constrained to be a power of two and therefore its objective function value is

[TABLE]

Consider now a scenario where each scaling transform was initialized and each $\mathbf{S}_{i_{k}}$ is updated again to further reduce the objective function while all others are kept fixed. Similarly to (21), the objective function develops now to

[TABLE]

where $\mathbf{A}_{k}=\prod_{q=k+1}^{m}\mathbf{S}_{i_{q}}$, $\mathbf{X}_{k}=\prod_{q=1}^{k-1}\mathbf{S}_{i_{q}}\mathbf{X}$ and the structure of $\mathbf{S}_{i_{k}}=\mathbf{I}_{n\times n}+(a_{k}-1)\text{diag}(\mathbf{e}_{i_{k}})$. For the other variables we have used here the same notation as in (21). Since we want to minimize the quantity in (29), then we have

[TABLE]

and using the minimizer of this expression

[TABLE]

it follows that for $i=1,\dots,n$ we have

[TABLE]

where $\mathbf{b}_{i_{k}}=\mathbf{x}_{i_{k}}\otimes\mathbf{a}_{i_{k}}$ is the $i_{k}^{\text{th}}$ column of $\mathbf{B}_{k}$ and $\mathbf{x}_{i_{k}}^{T}$ is the $i_{k}^{\text{th}}$ row of $\mathbf{X}_{k}$ and $\mathbf{a}_{i_{k}}$ is the $i_{k}^{\text{th}}$ column of $\mathbf{A}_{k}$. We have used the fact that the squared $\ell_{2}$ norm of a Kronecker product is the product of the squared $\ell_{2}$ norms of the two vectors involved.

Due to their simplicity, the scaling transformations $\mathbf{S}_{i}$ are numerically efficient, i.e., matrix-vector multiplication $\mathbf{S}_{i}\mathbf{x}$ takes only one bit shift or one multiplication operation (if we set $\alpha^{\star}=1$ and just solve a general unconstrained least squares problem on the $i^{\text{th}}$ coordinate). The inverse $\mathbf{S}_{i}^{-1}$ is also a scaling transformation (with $a^{-1}$ on position $i$) with the same numerical properties as the direct transformation.

To encode one scaling transformation we need $C+\log_{2}n$ bits (the constant is the cost of encoding of $a$, say $C=64$ for a double float, while the second term encodes the index $i$). Notice that we could represent $a$ in $\mathcal{R}_{p}$, for a given $p$, but then $a^{-1}$ does not have, in general, a representation in $\mathcal{R}_{p^{\prime}}$ for any finite $p^{\prime}$. It is for this reason that when it comes to the scaling transformation and its scaling factor $a$ we allow only two possibilities: either we take $a^{\star}$ the optimum least squares choice if we are working with arbitrary precision, i.e., $p=\infty$, or we take $\alpha^{\star}a^{\star}$ to ensure the scaling is a bit shift and therefore ensure the consistency of the inverse scaling transformation.

V-C Numerically efficient general transforms: Sm,p–DLA

Similarly to the other transformations described in this paper, we consider now the following dictionary structure

[TABLE]

where each $\mathbf{S}_{i_{k}j_{k}}$ is either a scaling or a shear transformation. For convenience, we denote $m_{1}$ and $m_{2}$ the number of scalings and shears, respectively – we have $m_{1}+m_{2}=m$.

The complete learning method is described in Algorithm 4. This procedure has two main components: the initialization phase and the iterative process that improves the factorization.

In the initialization phase, we use again the singular value decomposition, assume that our initial dictionary is $\mathbf{U}$ and then proceed to compute the first sparse representations $\mathbf{X}$. Then we proceed to initialize each $\mathbf{S}_{i_{k}j_{k}}$ iteratively. The computational cost of this step is dominated by the calculation of all the scores which takes $O(n^{3}N)$: there are $n^{2}$ scores and we need $O(nN)$ to compute each one (it is the computational cost of performing the dot product between $\mathbf{f}_{k}$ and the Kronecker product). Then, because we are dealing with a general dictionary, we use the batch Orthogonal Matching Pursuit (OMP) algorithm [31, 32] to build $\mathbf{X}$.

The iterative process now tries to improve the accuracy of the factorization by updating each individual transformation while all others are kept fixed. When we update a transformation we calculate the indices $(i_{k},j_{k})$ and the coefficients of the transformation $a,b$ or $c$. For this reason, this iterative step is computationally expensive at takes $O(n^{3}mNK)$: there are $K$ iterations which update each of the $m$ transformations and for each we need to compute $n^{2}$ scores where the computing load for a single one is dominated a dot product between $\mathbf{f}_{k}$ and a Kronecker product which takes $O(nN)$. Luckily, in practice we observer that a low number of iterations will usually suffice to reach a very low representation error $\|\mathbf{Y}-\mathbf{SX}\|_{F}^{2}$. After all the transformations are updated, we also update the representations $\mathbf{X}$ again by the OMP algorithm.

We allow an input parameter $p$ to set the precision of the transformations we compute. In the case of the shear transforms, we will represent the coefficients $b$ and $c$ in the set $\mathcal{R}_{p}$ while for the scaling transformations we allow two options: either use the full precision $a$ computed for $p=\infty$ or approximate it by the nearest power of two whenever $p\neq\infty$. We treat the scaling transformation in this binary fashion in order to keep consistent the inverse scaling operation (see the discussion in Section V-B).

Iteratively, we update each $\mathbf{S}_{i_{k}j_{k}}$ component of $\mathbf{S}$ and the sparse representations $\mathbf{X}$ via the OMP algorithm. Although the optimization steps to update the dictionary components are solved exactly and therefore always reduce the objective function, the OMP step cannot be guaranteed to reach optimal representations in general and therefore the algorithm is generally not guaranteed to converge monotonically to a local minimum point. As such, we track the best solution obtained so far in the iterative process and return it.

Transforms built by Sm,∞–DLA will perform $m$ multiplications and $m_{2}$ additions and those built by Sm,p–DLA will perform $m_{1}+pm_{2}$ bit shifts and $pm_{2}$ additions when $p\neq\infty$. Therefore, the computational complexity of these transforms is not predetermined only by the choice of $m$ but also by the actual factorization which is constructed. Related to our previously stated desired computational properties we have that: P1 is achieved by taking $m\sim O(\log n)$ and since inverses of scalings and shears are themselves scalings and shears (this covers also P3) we have the same computational benefits for the inverse transformation; for P2 and P4 we need the scalings to perform shift operations and the shears to have $p<\infty$.

Remark 6 (Completeness of scalings and shears). Every invertible $\mathbf{S}$ can be represented as a product of $n^{2}-n$ shear and $n$ scaling transformations (note that this was not the case for the binary orthogonal building blocks, see Remark 2). The proof of this fact is constructive: perform Gaussian elimination on $\mathbf{S}$ to diagonalize it using $n^{2}-n$ shears and then represent the resulting diagonal with $n$ scalings (permutations are also allowed if elimination is done with pivoting). Of course, in our approach we are trying to build $\mathbf{S}$ such that $m\ll n^{2}$, otherwise there is no computational benefit. $\hfill\blacksquare$

V-D A note on the general overcomplete case

Throughout this paper we have considered learning square linear transformations. Many of the successful applications of dictionary learning involves overcomplete transformations, i.e., learning an $n\times d$ transformation with $d>n$. All building blocks that we have introduced are naturally square, e.g., the factorization (33), but they can be made overcomplete by pre-multiplication with a mask matrix $\mathcal{M}=\begin{bmatrix}\mathbf{I}_{n\times n}&\mathbf{0}_{n\times(d-n)}\end{bmatrix}$. From an optimization perspective we are therefore in the situation of the quantities (21) and (29) where we replace $\mathbf{A}_{k}$ by $\mathcal{M}\mathbf{A}_{k}$. The approach to improve each building block is essentially the same with the notable exception that many scores are null.

VI Experimental results

In this section we test the proposed learning algorithms on image data, where we have well-known numerically efficient transformations against which to compare. The training data we consider is built from popular test images from the image processing literature (lena, peppers, boat etc.). The dataset $\mathbf{Y}\in\mathbb{R}^{64\times 12288}$ consists of $8\times 8$ non-overlapping image patches with their means removed. To evaluate the learning algorithms, we take the relative representation error of the dataset $\mathbf{Y}$ in the dictionary $\mathbf{D}$ with the representations $\mathbf{X}$ as

[TABLE]

The dictionary learning problem constraints require that the transformation we learn have $\ell_{2}$ normalized columns. This constraint is trivially obeyed for the orthonormal transforms designed via Bm–DLA and Mq–DLA but not in the case of those built by Sm,p–DLA. One solution is to introduce a diagonal matrix $\mathbf{D}\in\mathbb{R}^{n\times n}$ and update (15) via $\mathbf{S}\leftarrow\mathbf{SD}$ such that columns are $\ell_{2}$ normalized. Once the representations are computed, $\mathbf{S}$ is restored to (33) and $\mathbf{X}\leftarrow\mathbf{DX}$.

We compare against the FFT or the discrete cosine transform (DCT) which, for real-valued inputs, has complexity $2n\log_{2}n$: $\frac{3}{2}n\log_{2}n$ additions and $\frac{1}{2}n\log_{2}n$ multiplications. We also compare against transforms built by Gm–DLA, which have complexity $6m$, even additions and multiplications.

In Figure 1 we show the representation errors achieved by Bm–DLA and compare them against Gm–DLA [1]. We always expect Gm–DLA to perform better, as shown in the figure, especially for large $m$. Since up to $m=64$ we have that Bm–DLA closely tracks Gm–DLA, this introduces the idea of a potential hybrid algorithm: use binary transformations in the beginning and then proceed with Gm–DLA when the decrease in the error slows or plateaus. The advantage of Bm–DLA transforms is that they are $33\%$ faster than Gm–DLA transforms. Similar results are shown in Figure 2 for Mq–DLA. The same plateau of the error is observed for $m\geq 128$ for all sparsity levels $s$. Because this approach is computationally expensive in the training phase, we also show a cheaper method called Mq–DLA (see Remark 3) that performs similarly. We complete the orthonormal learning experiments with Figure 3 where we show B${}_{m}^{\otimes}$–DLA (see Remark 4) for $s\in\{4,8,12\}$. Notice that the representation error drops faster than that of Bm–DLA and it plateaus for higher $m$. As explained, the disadvantage of B${}_{m}^{\otimes}$–DLA is the increased training time: on a modern computer, calculating all scores $C_{i_{k}j_{k}r_{k}p_{k}}^{(t_{k})}$ takes over three hours and updating the scores with each iteration is done only for a random subset to keep the running time for one iteration to only a few minutes.

In Figure 4 we show the representation errors achieved by Sm,p–DLA. The effect of the precision parameter $p$ is interesting. For any sparsity level $s$, when the number of transformations $m$ is low ($\leq 64$) the parameter $p$ does not play an important role and therefore $p=1$, which has the lowest computational complexity, is preferred. For other $m$, there are slight differences (mostly within $5\%$) in the error and, as expected, $p=\infty$ works best. For $s=4$ the proportion of scalings and shears in the transforms designed is approximately $10\%$ to $90\%$.

Remark 7 (Sub-optimal transform learning). All proposed algorithms choose indices $i,j$ on which the linear transformation operate greedily (maximally reduce the current objective function). In some situations, it might be convenient to make sub-optimal choices that take into account other goals.

For example, from a computational perspective, we might consider a highly local computational model, i.e., our algorithms perform operations, in place, on memory locations that are close such that they can exploit benefits of hierarchically memory structures (cache-oblivious algorithms [33]). In our case this might translate in constraints as $|i_{k}-j_{k}|\leq\epsilon$, i.e., for a particular transformation control the distance between the operands, and $|i_{k}-i_{k+1}|\leq\epsilon$ (and similarly for $j_{k}$, $j_{k+1}$), i.e., consecutive transformations operate on neighboring regions of memory. Another example is bounding the dynamic range [14, Section IV] of the intermediate stages of our transforms. For example, in the case of Bm–DLA, given an input $\mathbf{x}$ such that $|x_{i}|\leq C$ the output $\mathbf{y}=\mathbf{Bx}$ is such that $|y_{i}|\leq 2^{\frac{m}{2}}C$ (equality happens if the same operation takes places for the same indices $i,j$, e.g., $x_{i}\leftarrow 2^{-\frac{1}{2}}(x_{i}+x_{j})$ happens $m$ times. Therefore, imposing some index diversity and/or making sure we do not have repetitive operations for the same indices will reduce the maximum possible output magnitude.$\hfill\blacksquare$

Finally, in Figures 5 and 6 we show a comparative study of the proposed transforms for fixed $s=4$. The point is to have a Pareto curve to show the computational-representation error trade-off: in the first plot we show the overall computational complexity and then, in the second plot, we apply an extra penalty $\gamma=6$ to multiplication operations. Unsurprisingly, the non-orthonormal transforms built by Sm,p–DLA are the most effective (for $p=\infty$ in Figure 5 and $p=1$ in Figure 6). The orthonormal transforms behave as expected as well: Mq–DLA performs very well in terms of computational complexity but worse in terms of representation errors (and plateaus quickly), Bm–DLA makes better progress and achieves better representation errors (although it also slows down progress with increased $m$) while B${}_{m}^{\otimes}$–DLA combines the benefits of the two approaches with the drawback is that training took overnight. Note that, for the same complexity, Bm–DLA and B${}_{m}^{\otimes}$–DLA are very close to the DCT while with the extra penalty $\gamma$ for multiplications B${}_{m}^{\otimes}$–DLA outperforms the DCT, i.e., when addition and bit shift operations are cheap enough we can afford a large enough number of these operations to surmount the constraint of avoiding multiplications.

VII Conclusions

In this paper we have proposed several dictionary learning algorithms that produce linear transformation which exhibit low computational complexity in general and reduce (or completely eliminate) multiplications in particular. We show that these transforms perform very well on image data where we compare against the DCT. We show that multiplications can be avoided and representation errors of image data can be kept low by replacing multiplication operations with a larger number of addition and bit shift operations.