Dictionary Learning with BLOTLESS Update

TL;DR

This paper introduces BLOTLESS, a novel block total least squares algorithm for dictionary update in sparse representation, which improves recovery accuracy especially with limited training data.

Contribution

The paper proposes BLOTLESS, a new dictionary update method that updates blocks of elements simultaneously, with theoretical recovery conditions and advantages over existing algorithms.

Findings

BLOTLESS achieves accurate dictionary recovery with fewer training samples.

Theoretical bounds closely match empirical data requirements.

Numerical experiments show improved performance over state-of-the-art methods.

Abstract

Algorithms for learning a dictionary to sparsely represent a given dataset typically alternate between sparse coding and dictionary update stages. Methods for dictionary update aim to minimise expansion error by updating dictionary vectors and expansion coefficients given patterns of non-zero coefficients obtained in the sparse coding stage. We propose a block total least squares (BLOTLESS) algorithm for dictionary update. BLOTLESS updates a block of dictionary elements and the corresponding sparse coefficients simultaneously. In the error free case, three necessary conditions for exact recovery are identified. Lower bounds on the number of training data are established so that the necessary conditions hold with high probability. Numerical simulations show that the bounds approximate well the number of training data needed for exact dictionary recovery. Numerical experiments further demonstrate several benefits of dictionary learning with BLOTLESS update compared with state-of-the-art algorithms especially when the amount of training data is small.

Tables

Table 1. TABLE I: Runtime (seconds per iteration of the dictionary learning process) comparison for different total least squares methods: m = l = 64 𝑚 𝑙 64 m=l=64 , θ = 5 / 64 𝜃 5 64 \theta=5/64 , and number of iterations n it = 50 subscript 𝑛 it 50 n_{\rm it}=50 . Results are averages of 50 independent trials using Matlab 2018b on a MacBook Pro with 16GB RAM and 2.3 GHz Intel Core i5 processor.

	$n=200$	$n=300$	$n=400$	$n=500$
STLS	621.6544	836.3027	1098.5955	1265.1895
ParTLS	9.2544	18.1191	25.3678	31.0551
IterTLS	6.3446	8.5157	10.9026	13.8056

Table 2. TABLE II: Denoising comparison using different learnt dictionaries, where the denoised PSNR (dB) are computed and shown in table.

Original Image
Noisy Image	28.13 dB	22.11 dB	18.59 dB	28.13 dB	22.11 dB	18.59 dB	28.13 dB	22.11 dB	18.59 dB	28.13 dB	22.11 dB	18.59 dB
MOD	33.00 dB	29.24 dB	27.22 dB	32.68 dB	28.84 dB	26.82 dB	31.95 dB	27.43 dB	25.76 dB	32.13 dB	28.23 dB	26.26 dB
K-SVD	32.50 dB	28.72 dB	26.79 dB	32.03 dB	28.26 dB	26.35 dB	31.49 dB	27.43 dB	25.31 dB	31.58 dB	27.74 dB	25.81 dB
SimCO	33.43 dB	29.78 dB	27.81 dB	33.58 dB	30.11 dB	28.04 dB	32.18 dB	28.51 dB	26.47 dB	32.65 dB	29.18 dB	27.27 dB
BLOTLESS	33.67 dB	29.90 dB	27.95 dB	33.91 dB	30.33 dB	28.20 dB	32.38 dB	28.67 dB	26.66 dB	32.88 dB	29.38 dB	27.45 dB

Equations

$$\underset{\mathbf{D},\mathbf{X}}{\min }\Vert \mathbf{Y}-\mathbf{D}\mathbf{X}{\Vert }_F^2,s.t.\Vert {\mathbf{X}}_{:,j}{\Vert }_0\le k,\forall j\in [n],$$

$$\underset{{\mathbf{X}}_{:,j}}{\min }\Vert {\mathbf{Y}}_{:,j}-\mathbf{D}{\mathbf{X}}_{:,j}{\Vert }_2,s.t.\Vert {\mathbf{X}}_{:,j}{\Vert }_0\le k,\forall j\in [n],$$

$$\underset{\mathbf{D}}{\min }\Vert \mathbf{Y}-\mathbf{D}\mathbf{X}{\Vert }_F^2.$$

$$\underset{\mathbf{D},\mathbf{X}}{\min }\Vert \mathbf{Y}-\mathbf{D}\mathbf{X}{\Vert }_F^2,\text{s.t.}{\mathcal{P}}_{{\Omega }^c}(\mathbf{X})=0$$

$$\underset{{\mathbf{D}}_{:,i},{\mathbf{X}}_{i,:}}{\min }\Vert {\left(\mathbf{Y}-{\mathbf{D}}_{:,\{i{\}}^c}{\mathbf{X}}_{\{i{\}}^c,:}\right)}_{:,{\Omega }_i}-{\mathbf{D}}_{:,i}{\mathbf{X}}_{i,{\Omega }_i}{\Vert }_F^2.$$

$$\underset{\mathbf{D}}{\min }\underset{f(\mathbf{D})}{\underbrace{\underset{\mathbf{X}:{\mathcal{P}}_{{\Omega }^c}(\mathbf{X})=0}{\min }\Vert \mathbf{Y}-\mathbf{D}\mathbf{X}{\Vert }_F^2}}.$$

$$\underset{{\mathbf{D}}_{:,\mathbf{i}}}{\min }\underset{{\mathbf{X}}_{\mathbf{i},:}:{\mathcal{P}}_{{\Omega }_i^c}({\mathbf{X}}_{i,:})=0}{\min }\Vert \mathbf{Y}-\mathbf{D}\mathbf{X}{\Vert }_F^2.$$

$$\mathbf{Y}=\mathbf{D}\mathbf{X}\mathrm{and}{\mathcal{P}}_{{\mathrm{\Omega }}^{\mathrm{c}}}(\mathbf{X})=0.$$

$$\mathbf{H}\mathbf{Y}=\mathbf{X}\mathrm{and}{\mathcal{P}}_{{\mathrm{\Omega }}^{\mathrm{c}}}(\mathbf{X})=0,$$

$$\left[\bm{Y}^{T},-\bm{I}_{n\times n}\right]\left[\begin{array}[]{c}\bm{H}^{T}\\ \bm{X}^{T}\end{array}\right]=\bm{0}_{n\times m}\;\;{\rm and}\;\;\mathcal{P}_{\Omega^{c}}(\bm{X})=\bm{0},$$

$$\mathbf{H}[\mathbf{Y},1_{m\times 1}]=[\mathbf{X},1_{m\times 1}]\mathrm{and}{\mathcal{P}}_{{\mathrm{\Omega }}^{\mathrm{c}}}(\mathbf{X})=0.$$

$$({\mathbf{H}}_0{)}_{i,:}\mathbf{Y}=({\mathbf{X}}_0{)}_{i,:},$$

$${\mathbf{H}}_{i,:}{\mathbf{Y}}_{:,{\Omega }_i^c}=0^T,\mathrm{and}{\mathbf{H}}_{i,:}{\mathbf{Y}}_{:,j}\ne 0\forall j\in {\Omega }_i.$$

$$n_1$$

$$+\sqrt{{\left(1-\frac{\ln ϵ}{4m\left(m-1\right)\left(1-\theta \right)}\right)}^2-1}].$$

$$n_2$$

$$+\sqrt{{\left(1-\frac{\ln ϵ-\ln m}{4\left(m-1\right)\left(1-\theta \right)}\right)}^2-1}].$$

$$n_3=\frac{\ln ϵ-\ln m-\ln \left(m-1\right)}{\ln \left(1-\theta \left(1-\theta \right)\right)}.$$

$$\mathbf{H}[\mathbf{Y},1]\approx [\mathbf{X},1]\mathrm{and}{\mathcal{P}}_{{\Omega }^c}(\mathbf{X})\approx 0.$$

$$\underset{\tilde{\mathbf{A}},\tilde{\mathbf{B}},\mathbf{X}}{\min }\Vert [\mathbf{A}-\tilde{\mathbf{A}},\mathbf{B}-\tilde{\mathbf{B}}]{\Vert }_F^2,\mathrm{subject}\mathrm{to}\tilde{\mathbf{A}}\mathbf{X}=\tilde{\mathbf{B}}.$$

$$\underset{\tilde{\mathbf{Y}},\tilde{\mathbf{X}},\mathbf{H}}{\min }$$

$$s.t.$$

$$\mathcal{L}(\tilde{\mathbf{Y}},\tilde{\mathbf{X}},\mathbf{H}):=\mathbf{H}[\tilde{\mathbf{Y}},1]-[\tilde{\mathbf{X}},1].$$

$$\mathcal{L}(\tilde{\mathbf{Y}},\tilde{\mathbf{X}},\mathbf{H})=\mathcal{L}(\hat{\mathbf{Y}},\hat{\mathbf{X}},\hat{\mathbf{H}})+\mathbf{J}(\mathbf{z}-\hat{\mathbf{z}}),$$

$$\underset{\tilde{\mathbf{Y}},\tilde{\mathbf{X}},\mathbf{H}}{\min }$$

$$s.t.$$

$${\mathbf{H}}_{i,:}[\mathbf{Y},1]\approx [{\mathbf{X}}_{i,:},1]\mathrm{and}{\mathbf{X}}_{i,{\Omega }_i^c}\approx 0,\forall i\in [l].$$

$$\underbrace{\left[\begin{array}[]{cc}\bm{Y}^{T}&-\bm{P}_{\Omega_{i}}\\ \bm{1}&\bm{0}\end{array}\right]}_{\bm{A}_{i}}\underbrace{\left[\begin{array}[]{c}\bm{H}_{i,:}^{T}\\ \bm{X}_{i,\Omega_{i}}^{T}\end{array}\right]}_{\bm{z}_{i}}\approx\underbrace{\left[\begin{array}[]{c}\bm{0}\\ 1\end{array}\right]}_{\bm{b}},\;\forall i\in\left[l\right],$$

$$\underset{\tilde{\mathbf{Y}},\tilde{\mathbf{X}},\mathbf{H}}{\min }{\Vert [{\mathbf{Y}}^T-{\tilde{\mathbf{Y}}}^T,{\hat{\mathbf{X}}}^T-{\tilde{\mathbf{X}}}^T]\Vert }_F^{2}s.t.{\tilde{\mathbf{Y}}}^T{\mathbf{H}}^T={\tilde{\mathbf{X}}}^T,$$

$$\left[\bm{Y}^{T},\hat{\bm{X}}^{T}\right]=\left[\bm{U}_{Y},\bm{U}_{X}\right]\left[\begin{array}[]{cc}\bm{\Sigma}_{Y}&\bm{0}\\ \bm{0}&\bm{\Sigma}_{X}\end{array}\right]\left[\begin{array}[]{cc}\bm{V}_{YY}&\bm{V}_{YX}\\ \bm{V}_{XY}&\bm{V}_{XX}\end{array}\right]^{T}.$$

$${\mathbf{Y}}_r=\mathbf{Y}-{\mathbf{D}}_{:,{\mathcal{T}}^c}{\mathbf{X}}_{{\mathcal{T}}^c,:},$$

$$R_{\mathrm{err}}=\frac{1}{l}\sum _{p=1}^l\left(1-|{\hat{\mathbf{d}}}_p^T{\mathbf{d}}_{0_{j_p}}|\right),$$

Full text

I Introduction

Sparse signal representation has found a wide range of applications, including image denoising [1, 2], image in-painting[1], image deconvolution [3], image super-resolution [4, 5], etc. The key idea behind the concept of sparse representation is that natural signals tend to have sparse representations under certain bases/dictionaries. Hence, finding a dictionary under which a given data set can be represented in a sparse manner, has become a very active area of research. Although numerous analytical dictionaries exist, including Fourier basis[6], discrete cosine transform (DCT) dictionaries, wavelets[7], curvelets [8], etc., the need to adapt to properties of specific data sets has long been driving research efforts towards the development of efficient algorithms for dictionary learning [9, 10]. More formally, dictionary learning is the problem of finding a dictionary $\bm{D}\in\mathbb{R}^{m\times l}$ of $l$ vectors in $\mathbb{R}^{m}$ such that the $n$ training samples in $\bm{Y}\in\mathbb{R}^{m\times n}$ can be written as $\bm{Y}=\bm{D}\bm{X}$, where the coefficient matrix $\bm{X}\in\mathbb{R}^{l\times n}$ is sparse. Of particular interest is overcomplete dictionary learning where the number of dictionary items is larger than the data dimension, i.e., $l>m$, and the number of the training samples is typically much larger than the size of the dictionary, $n\gg l$. Dictionary learning is a nonconvex bilinear inverse problem, very challenging to solve in general.

The bilinear dictionary learning problem is typically approached by alternating between two stages: sparse coding and dictionary update [9, 10, 11, 12, 13, 14]. In the sparse coding stage, the goal is to find sparse representations $\bm{X}$ of training samples $\bm{Y}$ for a given dictionary $\bm{D}$. For that purpose, scores of algorithms have been developed. They can be divided into two main categories. The first category consists of greedy algorithms, including orthogonal matching pursuit (OMP) [15], regularized orthogonal matching pursuit (ROMP) [16], subspace pursuit (SP) [17], etc. In the second category, sparse coding is formulated as a convex optimization problem where $\ell_{1}$-norm is used to promote sparsity [18], and then optimization techniques, e.g. the fast iterative shrinkage-thresholding algorithm (FISTA) [19], can be applied. Reviews of sparse recovery algorithms can be found in [20].

The goal of the dictionary update is to refine the dictionary so that the training samples $\bm{Y}$ have more accurate sparse representations given indices of non-zero coefficients obtained in the sparse coding stage. In the probabilistic framework, one may apply either maximum likelihood (ML) estimator [9] or maximum a posteriori (MAP) estimator [12], and then solve them by using gradient decent procedures. In the context of ML formulation [9], Engan et al.  [11] proposed the method of optimal directions (MOD) where the sparse coefficients $\bm{X}$ are fixed and the dictionary update problem is cast as a least squares problem which can be solved efficiently; modifications of MOD were subsequently proposed in [21, 22, 23].

Recently an alternative approach for dictionary update has become popular, where both the dictionary and the sparse coefficients are updated simultaneously with a given sparsity pattern. The representative algorithms include the famous K-SVD algorithm [10, 24] and SimCO [13]. The crux of K-SVD [10] algorithm is to update dictionary items and their corresponding sparse coefficients simultaneously, sequentially one by one. K-SVD was subsequently extended to allow simultaneous update of multiple dictionary elements and corresponding coefficients [24]. SimCO [13], of which K-SVD is a special case, goes further and updates the whole dictionary and sparse coefficients simultaneously. The main idea of SimCO is that given a sparsity pattern, the sparse coefficients can be viewed as a function of the dictionary. As a result, the dictionary update becomes a nonconvex optimisation problem with respect to the dictionary. The optimisation is then preformed using the gradient descent method combined with a heuristic sub-routine designed to deal with singular points which can prevent from the convergence to the global minimum[13].

Due to the non-convexity of dictionary learning problem, it is challenging to understand under which conditions exact dictionary recovery is possible and which method is optimal in achieving that. Following early efforts on theoretical analysis of exact dictionary recovery [25, 26, 27, 28, 29, 30, 31], more recently, Spielman et. al. [32] studied dictionary learning problem with complete dictionaries where the dictionary can be presented as a square matrix. By solving a certain sequence of linear programs, they showed that one can recover a complete dictionary $\bm{D}$ from $\bm{Y}=\bm{D}\bm{X}$ when $\bm{X}$ is a sparse random matrix with $O(\sqrt{m})$ nonzeros per column. In [33, 34, 35, 36], the authors propose algorithms which combine clustering, spectral initialization, and local refinement to recover overcomplete and incoherent dictionaries.

Again these algorithms succeed when $\bm{X}$ has $O(\sqrt{m})$ nonzeros per column. The work in [37] provides a polynomial-time algorithm that recovers a large class of over-complete dictionaries when $\bm{X}$ has $O(m^{1-\delta})$ nonzeros per column for any constant $\delta\in(0,1)$. However, the proposed algorithm runs in super-polynomial time when the sparsity level goes up to $O(m)$. Similarly, in [38] the authors proposed a super-polynomial time algorithm that guarantees recovery with close to $O(m)$ nonzeros per column. Sun et al. [39, 40], on the other hand, proposed a polynomial-time algorithm that provably recovers complete dictionary $\bm{D}$ when $\bm{X}$ has $O(m)$ nonzeros per column and the size of training samples is $O(m^{2}\log(m))$.

This paper addresses the dictionary update problem, where both the dictionary and the sparse coefficients are updated, for a given sparsity pattern. Whilst it is a sub-problem of the overall dictionary learning, it is nevertheless an important step towards solving the overall problem, and its bilinear nature makes it nonconvex and hence very challenging to solve. Our main contributions are as follows.

• •

BLOTLESS simultaneously updates a block of dictionary items and the corresponding sparse coefficients. Inspired by ideas presented in [41, 42], in BLOTLESS the bilinear nonconvex block update problem is transformed into a linear least squares problem, which can be solved efficiently.

• •

For the error-free case, when the sparsity pattern is known exactly, three necessary conditions for unique recovery are identified, expressed in terms of lower bounds on the number of training data. Numerical simulations show that the theoretical bounds well approximate the empirical number of training data needed for exact dictionary recovery. In particular, we show that the number of training samples needed is $O(m)$ for complete dictionary update.

• •

BLOTLESS is numerically demonstrated robust to errors in the assumed sparsity pattern. When embedded into the overall dictionary learning process, it leads to faster convergence rate and less training samples needed compared to state-of-the-art algorithms including MOD, K-SVD and SimCO.

Our work is inspired by a recent work [41] where bilinear inverse problems are formulated as linear inverse problems. The main difference is that our theoretical analysis and algorithm designs in Sections III and IV are specifically tailored to the generic dictionary update problem while the focus of [41] is self-calibration which can be viewed as dictionary learning with only diagonal dictionaries. Parts of the results in this paper were presented in the conference paper [43]. In this journal paper, we refine the bounds in Section III and provide detailed proofs, add two total least squares algorithms in Section IV, and include more simulation results in Section V to support the design of the algorithm.

This paper is organized as follows. Section II briefly reviews dictionary learning and update methods. Section III discusses an ideal case where exact dictionary recovery is possible, for which a least squares method is developed and analysed. In Section IV, the general case of dictionary update is discussed, and the least squares method is extended to total least squares methods, leading to BLOTLESS. Results of extensive simulations are presented in Section V and conclusions are drawn in Section VI.

I-A Notation

In this paper, $\|\cdot\|_{2}$ denotes the $\ell_{2}$ norm and $\|\cdot\|_{F}$ stands for the Frobenius norm. For a positive integer $n$, define $[n]=\{1,2,\cdots,n\}$. For a matrix $\bm{M}$, $\bm{M}_{i,:}$ and $\bm{M}_{:,j}$ denote the $i$-th row and the $j$-th column of $\bm{M}$ respectively. Consider the sparse coefficient matrix $\bm{X}$. Let $\Omega$ be the support set of $\bm{X}$, i.e., the index set that containing indices of all nonzero entries in $\bm{X}$. Let $\Omega_{i}$ be the support set of the row vector $\bm{X}_{i,:}$. Then $\bm{X}_{i,\Omega_{i}}$ is the row vector obtained by keeping the nonzero entries of $\bm{X}_{i,:}$ and removing all its zero entries. Symbols $\bm{I}$, $\bm{1}$, and $\bm{0}$ denote the identity matrix, the vector of which all the entries are 1, and the vector with all zero entries, respectively. For a given set $\Omega_{i}\subset[n]$, $\Omega_{i}^{c}$ denotes its complement in $[n]$.

II Dictionary Learning: The Background

Dictionary learning is the process of finding a dictionary which sparsely represents given training samples. Let $\bm{Y}\in\mathbb{R}^{m\times n}$ be the training sample matrix, where $m$ is the dimension of training sample vectors and $n$ is the number of training samples. The overall dictionary learning problem is often formulated as:

[TABLE]

where $\bm{D}\in\mathbb{R}^{m\times l}$ is the dictionary, $\bm{X}\in\mathbb{R}^{l\times n}$ is the sparse coefficient matrix, the $\ell_{0}$ pseudo-norm $\|\cdot\|_{0}$ gives the number of non-zero elements, also known as sparsity level, and $k<l$ is the upper bound of the sparsity level.

Dictionary learning algorithms typically iterate between two stages: sparse coding and dictionary update. The goal of sparsity coding is to find a sparse coefficient matrix $\bm{X}$ for a given dictionary $\bm{D}$. One way to achieve this is to solve the problem

[TABLE]

using greedy algorithms, for example OMP [15] or SP [17].

In the dictionary update stage, the goal is to refine the dictionary with either fixed sparse coefficients or a fixed sparse pattern, i.e. fixed locations of non-zero coefficients. The famous MOD method [11] falls into the first category. With fixed sparse coefficients, dictionary update is simply a least squares problem

[TABLE]

A more popular and advantageous approach is to simultaneously update the dictionary and nonzero sparse coefficients by fixing only the sparsity pattern. With this idea, dictionary update is then formulated as [10, 13, 24]

[TABLE]

where $\mathcal{P}_{\Omega^{c}}(\bm{X})$ gives the vector formed by the entries of $\bm{X}$ indexed by $\Omega^{c}$. However, problem (3) is bilinear, nonconvex, and challenging to solve.

Among many methods for solving (3), we here briefly review K-SVD [10] and SimCO [13]. K-SVD algorithm successively updates individual dictionary items $\bm{D}_{:,i}$ and the corresponding sparse coefficients $\bm{X}_{i,:}$ whilst keeping all other dictionary items and coefficients fixed:

[TABLE]

The optimal solution can be obtained by taking the largest left and right singular vectors of the matrix $\left(\bm{Y}-\bm{D}_{:,\{i\}^{c}}\bm{X}_{\{i\}^{c},:}\right)_{:,\Omega_{i}}$.

The idea of SimCO is to formulate the dictionary update problem in (3) as a nonconvex optimisation problem with respect to the overall dictionary, that is

[TABLE]

and then solve it using gradient descent of $\bm{D}$. This leads to an update of all dictionary items and sparse coefficients simultaneously. K-SVD can be viewed as a special case of SimCO where the objective function reads

[TABLE]

The focus of this paper is a novel solution to Problem (3).

III Exact Dictionary Recovery

This section focuses on an ideal case that the dictionary can be exactly recovered. We assume that the training samples in $\bm{Y}$ are generated from $\bm{Y}=\bm{D}_{0}\bm{X}_{0}$ where $\bm{D}_{0}$ is a tall or square matrix ($m\geq l$) and the sparsity pattern of $\bm{X}_{0}$ is given. For compositional convenience, we focus on the case where $\bm{D}_{0}$ is a square matrix, $\bm{D}_{0}\in\mathbb{R}^{m\times m}$, as the same analysis is valid for a tall dictionary where $m>l$.

With given sparsity pattern denoted by $\Omega$, the dictionary update problem can be formulated as a bilinear inverse problem in which the goal is to find $\bm{D}$ and $\bm{X}$ such that

[TABLE]

The constraint $\bm{Y}=\bm{D}\bm{X}$ is nonconvex. Generally speaking, it is challenging to solve (6) and there are no guarantees for the global optimality of the solution.

III-1 Least Squares Solver

Suppose that the unknown dictionary matrix $\bm{D}$ is invertible. The nonconvex problem in (6) can be translated into a convex problem by using a strategy similar to that explored in [41]. Define $\bm{H}=\bm{D}^{-1}$. Then $\bm{H}\bm{Y}=\bm{X}$. The goal is now to find $\bm{H}$ and $\bm{X}$ such that

[TABLE]

or equivalently,

[TABLE]

where the subscripts are used to indicate matrix dimensions. In this manner the original bilinear problem (6) is cast as an equivalent linear least squares problem.

However, the formulation in (8) admits trivial solution $\bm{H}=\bm{0}$ and $\bm{X}=\bm{0}$. In fact, (8) admits at least $m$ linearly independent solutions.

Proposition 1**.**

There are at least $m$ linear independent solutions to the least squares problem in (8).

Proof.

This proposition is proved by construction. Let $\bm{H}_{0}=\bm{D}_{0}^{-1}$. Define matrix $\bm{Z}_{i}\in\mathbb{R}^{\left(m+n\right)\times m}$ by keeping the $i$-th column of the matrix $\left[\bm{H}_{0},\bm{X}_{0}\right]^{T}$ and setting all other columns to zero, that is, $\left(\bm{Z}_{i}\right)_{:,i}=\left[(\bm{H}_{0})_{i,:},(\bm{X}_{0})_{i,:}\right]^{T}$ and $\left(\bm{Z}_{i}\right)_{:,j}=\bm{0}$ for all $j\neq i$. From the fact that $(\bm{H}_{0})_{i,:}\bm{Y}=(\bm{X}_{0})_{i,:}$, it is straightforward to verify that $\bm{Z}_{i}$, $i\in\left[m\right]$, is a solution of (8).

The solutions $\bm{Z}_{i}$, $i\in\left[m\right]$, are linearly independent. This can be easily verified by observing that the positions of nonzero elements in $\bm{Z}_{i}$ and $\bm{Z}_{j}$, $i\neq j$, are different. ∎

III-2 Necessary Conditions for Unique Recovery

We now consider the uniqueness of the solution in more detail and derive necessary conditions for unique recovery. Two ambiguities can be identified in the dictionary update problem in (8). The first is permutation ambiguity. Let $\Omega_{i}$ and $\Omega_{j}$ be the support sets (the index set containing indices corresponding to nonzero entries) of the $i$-th and $j$-th row of $\bm{X}_{0}$. If $\Omega_{i}=\Omega_{j}$, then the tuple $(\bm{D}_{0}\bm{P}_{i\Leftrightarrow j},\bm{P}_{i\Leftrightarrow j}\bm{X}_{0})$ is a valid solution of (6), where $\bm{P}_{i\Leftrightarrow j}$ denotes the permutation matrix generated by permuting the $i$-th and $j$-th row of the identity matrix. On the other hand, there is no permutation ambiguity if $\Omega_{i}\neq\Omega_{j}$ for all $i\neq j$. In practice, the given sparsity pattern is typically diverse enough to avoid permutation ambiguity.

The second is the scaling ambiguity which cannot be avoided. Let $\bm{S}$ be a diagonal matrix with nonzero diagonal elements. It is clear that the tuple $(\bm{D}_{0}\bm{S},\bm{S}^{-1}\bm{X}_{0})$ is also a valid solution of (6). All solutions of the form $(\bm{D}_{0}\bm{S},\bm{S}^{-1}\bm{X}_{0})$ form an equivalent class. The scaling ambiguity can be addressed by ntroducing additional constraints. One option used in [41] is that the sum of the elements in each row of $\bm{H}$ is one, i.e., $\sum_{j}\bm{H}_{i,j}=1,\forall i$. With these constraints, one has

[TABLE]

Henceforth, we define unique recovery as unique up to the scaling ambiguity.

Definition 1** (Unique Recovery).**

The dictionary update problem (6) is said to admit a unique solution if all solutions are of the form $\bm{D}=\bm{D}_{0}\bm{S}$ and $\bm{X}=\bm{S}^{-1}\bm{X}_{0}$ for some diagonal matrix $\bm{S}$ with nonzero diagonal elements.

In the following, we identified three necessary conditions for unique recovery.

Proposition 2**.**

Assume that $\bm{D}_{0}$ is square and invertible. If the problem (6) has unique solution, then it holds that

$n\geq n_{0}=m+\frac{|\Omega|}{m}-1$. 2. 2.

For all $i\in[m]$, the support set of the $i$-th row of $\bm{X}_{0}$, denoted by $\Omega_{i}$, satisfies $|\Omega_{i}^{c}|\geq m-1$. 3. 3.

For all $i\in[m]$ and all $i^{\prime}\neq i$, $\exists j\in\left[n\right]$ such that $(\bm{X}_{0})_{i,j}=0$ and $(\bm{X}_{0})_{i^{\prime},j}\neq 0$.

Proof.

Necessary condition 1 is proved by using the fact that the solution of (9) is unique only if the number of equations is larger or equal than the number of unknown variables. The number of unknown variables in (9) is $(n+m)m$ while the number of equations in (9) is $(n+1)m+(nm-|\Omega|)$. Elementary calculations lead to the bound $n_{0}$.

The proof of the other two necessary conditions is based on the fact that

[TABLE]

where $\bm{H}_{0}=\bm{D}_{0}^{-1}$. To simplify the notations, we omit the subscript 0 from $\bm{H}_{0}$, $\bm{D}_{0}$ and $\bm{X}_{0}$ in the rest of this proof. Divide the sample matrix $\bm{Y}$ into two sub-matrices $\bm{Y}_{:,\Omega_{i}}$ and $\bm{Y}_{:,\Omega_{i}^{c}}$. Then it holds that

[TABLE]

($\bm{H}_{i,:}^{T}$ is in the null space of $\bm{Y}_{:,\Omega_{i}^{c}}$.) Hence $\bm{H}_{i,:}$ is unique (up to a scaling factor) if and only if ${\rm span}(\bm{Y}_{:,\Omega_{i}^{c}})$ has dimension $m-1$. In this case, ${\rm span}(\bm{H}_{i,:}^{T})$ is the null space of both ${\rm span}(\bm{Y}_{:,\Omega_{i}^{c}})$ and ${\rm span}(\bm{D}_{:,\{i\}^{c}})$. It is concluded that $\bm{H}_{i,:}$ is unique if and only if ${\rm span}(\bm{Y}_{:,\Omega^{c}_{i}})={\rm span}(\bm{D}_{:,\left\{i\right\}^{c}})$.

Necessary condition 2 follows directly from that ${\rm rank}(\bm{Y}_{:,\Omega^{c}_{i}})=m-1$.

To prove the last necessary condition, note first that the fact that ${\rm span}(\bm{Y}_{:,\Omega^{c}_{i}})={\rm span}(\bm{D}_{:,\left\{i\right\}^{c}})$ implies that each column of $\bm{D}_{:,\left\{i\right\}^{c}}$ participates in generating some columns of $\bm{Y}_{:,\Omega^{c}_{i}}$. That is, $\forall i^{\prime}\neq i$, $\bm{D}_{:,i^{\prime}}$ participates in generating $\bm{Y}_{:,j}$ for some $j\notin\Omega^{c}_{i}$. Necessary condition 3 is therefore proved. Note that condition 3 is not sufficient. It does not prevent the following rank deficient case: there exist $i_{1}^{\prime},i_{2}^{\prime}\in\left\{i\right\}^{c}$ such that both $\bm{D}_{:,i_{1}^{\prime}}$ and $\bm{D}_{:,i_{2}^{\prime}}$ only participate in generating a single sample in $\bm{Y}_{:,j}$ for some $j\in\Omega_{i}^{c}$. ∎

III-3 Discussions on the Number of Samples

We now study the number of samples $n$ needed to ensure that the necessary conditions for unique recovery, as specified in Proposition 2, hold with high probability. To that end we use the following probabilistic model: entries of $\bm{D}_{0}\in\mathbb{R}^{m\times m}$ are independently generated from the Gaussian distribution $\mathcal{N}(0,\frac{1}{m})$, and entries of $\bm{X}_{0}\in\mathbb{R}^{m\times n}$ are independently generated from the Bernoulli-Gaussian distribution $BG(\theta)$ with $\theta\in[0,1]$, where Bernoulli-Gaussian distribution is defined as follows.

Definition 2**.**

A random variable $X$ is Bernoulli-Gaussian distributed $X\sim BG(\theta)$ with $\theta\in\left[0,1\right]$, if $X=W\cdot C$, where random variables $W$ and $C$ are independent, $W$ is Bernoulli distributed with parameter $\theta$, and $C\sim\mathcal{N}(0,1)$.

Remark 1*.*

The Gaussian distribution is not essential. It can be replaced by any continuous distribution.

Proposition 3** (Number of Samples).**

Suppose that $\bm{Y}=\bm{D}_{0}\bm{X}_{0}$ where $\bm{D}_{0}$ and $\bm{X}_{0}$ are generated according to the above probability model. Given a constant $\epsilon\in(0,1)$, the $i$-th necessary condition in Proposition 2 holds with probability at least $1-\epsilon$, if $n\geq n_{i}$, where

[TABLE]

[TABLE]

and

[TABLE]

Furthermore, it holds that $n_{1}\leq n_{2}$. If $n\geq\max(n_{2},n_{3})$, then all three necessary conditions in Proposition 2 hold.

Proof.

See Appendix A. ∎

Remark 2*.*

We have the following observations.

• •

With fixed $\epsilon$ and $\theta$, $n_{1}$ and $n_{2}$ scale linearly with $m$ while $n_{3}$ is proportional to $\ln m$.

• •

With fixed $m$ and $\theta$, $n_{1}$, $n_{2}$, and $n_{3}$ increase proportionally to $-\ln\epsilon$.

• •

With fixed $m$ and $\epsilon$, when $\theta$ increases from 0 to 1, $n_{1}$ and $n_{2}$ increase, while $n_{3}$ first decreases and then increases. This matches the intuition that when $\theta$ is too small, we need more samples to have enough information to recover the dictionary. On the other hand, when $\theta$ is too large, more samples are needed to generate the orthogonal space of each $\bm{H}^{T}_{i,:}$. This is verified by simulations in Section V.

The bound $\max(n_{2},n_{3})$ provides a good estimate of the number of samples needed for unique recovery. By set theory, if event $A$ is a necessary condition for $B$, then $B$ implies $A$, or equivalently, $B\subseteq A$ and $\Pr(B)\leq\Pr(A)$. In Proposition 3, the quantity $1-\epsilon$ is a lower bound for $\Pr(A)$, where these necessary conditions hold. But unfortunately it is neither lower nor upper bound for $\Pr(B)$, where the dictionary can be uniquely recovered. Nevertheless, our simulations show that $\max(n_{2},n_{3})$ is a good approximation to the number of samples needed to recover the dictionary uniquely with probability more than $1-\epsilon$.

In an asymptotic regime, the bounds can be simplified.

Corollary 1** (Asymptotic Bounds).**

Consider the same settings as in Proposition 3. For a given $\theta\in(0,1)$, let $m,n\rightarrow\infty$ with $\frac{n}{m}\rightarrow\bar{n}\in\mathbb{R}^{+}$. If $\bar{n}>\frac{1}{1-\theta}$, then all three necessary conditions in Proposition 2 holds with a probability arbitrary close to 1.

This corollary follows from elementary calculations and the fact that $\ln(m)/m\rightarrow 0$ when $m\rightarrow\infty$.

IV Dictionary Update with Uncertainty

While Section III studies the ideal case, this section focuses on the general case using the insight from Section III. In practice, there may be noise in the training samples $\bm{Y}$, and there may be errors in the assumed sparsity pattern. The exact equality in (6) may not hold any longer. Following the idea in Section III, total least squares methods are applied to handle the uncertainties. The techniques for non-overcomplete and overcomplete dictionaries are developed in Sections IV-A and IV-B respectively.

IV-A *Non-overcomplete Dictionary Update *

In the case $m\geq l$, let $\bm{H}=\bm{D}^{\dagger}$ be the pseudo-inverse of $\bm{D}$ and assume that $\bm{H}\bm{D}=\bm{I}_{l\times l}$. Due to the uncertainty, Equation (9) becomes approximate, that is,

[TABLE]

Total least squares is a technique to solve a least squares problem in the form $\bm{A}\bm{X}\approx\bm{B}$ where errors in both observations $\bm{B}$ and regression models $\bm{A}$ are considered[44, 45]. It targets at minimising the total errors via

[TABLE]

The constraint set above is nonconvex and hence (11) is a nonconvex optimisation problem. Nevertheless, its global optimal solution can be obtained by using the singular value decomposition (SVD). Set $\bm{Z}=[\bm{X}^{T},-\bm{I}]^{T}$. Observe that the constraint in (11) implies that $[\tilde{\bm{A}},\tilde{\bm{B}}]\bm{Z}=\bm{0}$. The optimal $\bm{Z}$ can be obtained from the smallest right singular vectors of the matrix $[\bm{A},\bm{B}]$, and the optimal $[\tilde{\bm{A}},\tilde{\bm{B}}]$ is a best lower-rank approximation of the matrix $[\bm{A},\bm{B}]$.

The difficulty in applying total least squares directly is due to the additional constraint $\mathcal{P}_{\Omega^{c}}(\mathbf{X})\approx\mathbf{0}$ in (10). Below we present three possible solutions, where the last one IterTLS excels and is adopted.

IV-A1 Structured Total Least Squares (STLS)

Consider having uncertainties in both $\bm{Y}$ and the sparsity pattern. Based on (10), a straightforward total least squares formulation is

[TABLE]

To solve the above nonconvex optimisation problem, we follow the approach in [46]. It involves an iterative process where each iteration solves an approximated quadratic optimisation problem which admits a closed-form optimal solution.

At each iteration, denote the initial estimate of $(\tilde{\bm{Y}},\tilde{\bm{X}},\bm{H})$ by $(\hat{\bm{Y}},\hat{\bm{X}},\hat{\bm{H}})$. Note that the constraint set in (12) can be written as $\mathcal{L}(\tilde{\bm{Y}},\tilde{\bm{X}},\bm{H})=\bm{0}$ where

[TABLE]

We consider the first order Taylor approximation of $\mathcal{L}(\tilde{\bm{Y}},\tilde{\bm{X}},\bm{H})$ at given $(\hat{\bm{Y}},\hat{\bm{X}},\hat{\bm{H}})$, which reads

[TABLE]

where $\bm{z}:=[{\rm vect}(\tilde{\bm{Y}})^{T},{\rm vect}(\tilde{\bm{X}})^{T},{\rm vect}(\bm{H})^{T}]^{T}$, $\hat{\bm{z}}:=[{\rm vect}(\hat{\bm{Y}})^{T},{\rm vect}(\hat{\bm{X}})^{T},{\rm vect}(\hat{\bm{H}})^{T}]^{T}$, and $\bm{J}$ is the corresponding Jacobian matrix. With this approximation, the nonconvex optimisation problem in (12) becomes a quadratic optimisation problem with equality constraints

[TABLE]

This is a quadratic optimisation problem with linear equality constraints, and admits a closed-form solution by a direct application of KKT conditions [47].

The STLS approach has two issues. The first issue is its very high computational cost. The quadratic optimisation problem (13) involves $mn+ln+lm$ unknowns and $l(n+1)$ equation constraints. Its closed-form solution involves a matrix of size $(n(m+2l)+l(m+1))\times(n(m+2l)+l(m+1))$. We have obtained the closed-form of the Jacobian matrix $\bm{J}$, implemented a conjugate gradient algorithm to use the structures in (13) for a speed-up (details are omitted here). However, simulations in Section V-B show that the computation speed is still too slow for practical problems. The second issue is the inferior performance compared to other TLS methods in Sections IV-A2 and IV-A3. This is because Taylor approximation of the constraint is used in STLS, while other TLS methods below incorporate the constraints directly without Taylor approximation.

IV-A2 Parallel Total Least Squares (ParTLS)

The key idea of ParTLS is to decouple the problem (10) into $l$ sub-problems that can be solved in parallel:

[TABLE]

It is straightforward to verify that this is equivalent to

[TABLE]

where $\bm{P}_{\Omega_{i}}\in\mathbb{R}^{n\times|\Omega_{i}|}$ is the projection matrix obtained by keeping the columns of the identity matrix indexed by $\Omega_{i}$ and removing all other columns.

Sub-problems (14) can be solved by directly applying the TLS formulation (11). Note that $[\tilde{\bm{A}}_{i},-\tilde{\bm{b}}][\bm{z}_{i}^{T},1]^{T}=\bm{0}$. The vector $[\bm{z}_{i}^{T},1]^{T}$ can be computed as a scaled version of the least right singular vector of the matrix $[\bm{A}_{i},-\bm{b}]$. Then $\bm{H}_{i,:}$ and $\bm{X}_{i,\Omega_{i}}$ can be obtained from $\bm{z}_{i}$.

ParTLS enjoys the following advantages. 1) Its global optimality is guaranteed for the ideal case of no data noise or sparsity pattern errors. It is straightforward to see that in the ideal case the ParTLS solutions satisfy (9). 2) It is computationally efficient. The sub-problems (14) are of small size and can be solved in parallel. However, ParTLS also has its own issue — certain structures in the problem are not enforced. For different sub-problem $i\in[l]$, the ‘denoised’ $\bm{Y}$, denoted by $\tilde{\bm{Y}}_{i}$ can be different.

IV-A3 Iterative Total Least Squares (IterTLS)

IterTLS is an iterative algorithm such that in each iteration a total least squares problem is formulated based on the estimate from the previous iteration. It starts with an initial estimate obtained by solving the ideal case equation (9). In each iteration, let $\hat{\bm{X}}$ be an estimate of $\bm{X}$ from either initialisation or the previous iteration. We formulate the following total least squares problem

[TABLE]

which has the identical form as (11). Note that the constraint $\mathcal{P}_{\Omega^{c}}(\bm{X})\approx\bm{0}$ in (10) is implicitly imposed as $\mathcal{P}_{\Omega^{c}}(\hat{\bm{X}})=\bm{0}$. The problem (15) can be optimally solved by using the SVD

[TABLE]

The optimal solution is given by $\tilde{\bm{Y}}^{T}=\bm{U}_{Y}\bm{\Sigma}_{Y}\bm{V}_{YY}^{T}$ and $\tilde{\bm{X}}^{T}=\bm{U}_{Y}\bm{\Sigma}_{Y}\bm{V}_{XY}^{T}$. To prepare the next iteration, one obtains an updated estimate $\hat{\bm{X}}^{\prime}$ by applying a simple projection operator to $\tilde{\bm{X}}$: $\mathcal{P}_{\Omega}(\hat{\bm{X}}^{\prime})=\mathcal{P}_{\Omega}(\tilde{\bm{X}})$ and $\mathcal{P}_{\Omega^{c}}(\hat{\bm{X}}^{\prime})=\bm{0}$. With this new estimate $\hat{\bm{X}}^{\prime}$, one can proceed with the next iteration until convergence.

IV-B *Update Overcomplete Dictionary *

The difficulty of overcomplete dictionary update comes from the fact that for an overcomplete dictionary $\bm{D}_{0}$, $\bm{D}_{0}^{\dagger}\bm{D}_{0}\neq\bm{I}$ in general, where $\bm{D}_{0}^{\dagger}$ is the pseudo-inverse of $\bm{D}_{0}$. Therefore, the above least squares or total least squares approaches cannot be directly applied.

To address this issue, a straightforward approach is to divide the whole dictionary into a set of sub-dictionaries each of which is either complete or undercomplete, and then update these sub-dictionaries one-by-one whilst fixing all other sub-dictionaries and the corresponding coefficients. More explicitly, given estimated $\bm{D}$ and $\bm{X}$, consider updating $\bm{D}_{:,\mathcal{T}}$, the submatrix of $\bm{D}$ consisting of columns indexed by $\mathcal{T}$, and $\bm{X}_{\mathcal{T,:}}$, the submatrix of $\bm{X}$ consisting of rows indexed by $\mathcal{T}$. Then, consider the residual matrix

[TABLE]

and apply the method in Section IV-A3 to solve the problem $\bm{Y}_{r}\approx\bm{D}_{:,\mathcal{T}}\bm{X}_{\mathcal{T},:}$. Then repeat this step for all sub-dictionaries. As the dictionary is updated block by block, we refer to our algorithm as BLOck Total LEast SquareS (BLOTLESS).

V Numerical Test

Parts of the numerical tests are based on synthetic data. The training samples $\bm{Y}=\bm{D}_{0}\bm{X}_{0}$ are generated according to the probability model specified in Section III-3. When the dictionary recovery is not exact, the performance criterion is the difference between the ground-truth dictionary $\bm{D}_{0}$ and the estimated dictionary $\hat{\bm{D}}$. In particular, the estimation error is defined as

[TABLE]

where $\bm{\hat{d}}_{p}$ is the $p$-th item in the estimated dictionary, $\bm{d}_{0_{j_{p}}}$ is the $j_{p}$-th item in the ground-truth dictionary, $j_{p}=\arg\max_{j\in\mathcal{J}_{p}}\{|\bm{\hat{d}}_{p}^{T}\bm{d}_{0_{j}}|\}$ and $\mathcal{J}_{p}=[l]\setminus\{j_{1},j_{2},\cdots,j_{p-1}\}$. The items in both dictionaries are normalised to have unit $\ell_{2}$-norm.

Numerical tests based on real data are presented in Section V-D2 for image denoising. The performance metric is the Peak Signal-to-Noise-Ratio (PSNR) of the denoised images.

V-A Simulations for Exact Dictionary Recovery

In this section we evaluate numerically bounds in Proposition 3. All the results presented here are based on 100 random and independent trials. For theoretical performance prediction, we compute $n_{2}$, $n_{3}$, and $\max(n_{2},n_{3})$ using $\epsilon=0.01$. In the numerical simulations, we vary $n$ and find its critical value $n_{\rm sim}$ under which exact recovery happens with an empirical probability at most 99% and above which exact recovery happens with an empirical probability at least 99%.

We start with the relation between the number of training samples $n$ and the sparsity ratio $\theta$ for a given $m$. In Figure 1, we fix $m=30$, vary $\theta$, and study the probability of exact recovery against the number of training samples. Results in Figure 1 show that the theoretical prediction $\max(n_{2},n_{3})$ is quite close to the critical $n_{\rm{sim}}$ obtained by simulations. One can also observe that the needed number of training samples first decreases and then increases as $\theta$ increases, which is also predicted by the theoretical bounds. A larger scale numerical test is done in Figure 2, where four sub-figures correspond to four different values of $m$. Once again, simulations demonstrate these bounds in Proposition 3 match the simulations very well.

Let us consider the required $n$ for exact recovery as a function of $m$ by fixing $\theta=0.2$. Simulation results are depicted in Figure 3. When $n_{3}\geq n_{2}$, $n_{\rm sim}$ behaves as $\ln m$. Otherwise, $n_{\rm sim}$ behaves as $m$. This is consistent with Proposition 3.

Finally, we numerically study the accuracy of the asymptotic results in Corollary 1. We draw normalised number of training samples $n/m$ for exact recovery as a function of sparsity ratio $\theta$ for different values of $m$, including $m=16,\;32,\;64$. Simulation results in Figure 4 show a trend that is consistent with the asymptotic results in Corollary 1.

V-B *Total Least Squares Methods *

The three total least squares methods introduced in Section IV are compared, henceforth denoted by BLOTLESS-STLS, BLOTLESS-ParTLS and BLOTLESS-IterTLS respectively, by being embedded in the dictionary learning process. Random dictionaries are used as the initial starting point of dictionary learning. OMP is used for sparse coding.

Fig. 5 compares the dictionary learning errors (17) for both complete and over-complete dictionaries. BLOTLESS-IterTLS converges the fastest and has the smallest error floor. Then a runtime comparison is given in Table I. BLOTLESS-IterTLS clearly outperforms the other two methods. It is therefore used as the default dictionary update method in later simulations.

V-C Robustness to Errors in Sparsity Pattern

Simulations are next designed to evaluate the robustness of different dictionary update algorithms to errors in sparsity pattern. Towards that end, a fraction $r$ of indices in the true support are randomly chosen to be replaced with the same number of randomly chosen indices not in the true support set. This erroneous sparsity pattern is then fed into different dictionary update algorithms. The numerical results in Fig. 6 demonstrate the robustness of BLOTLESS (with IterTLS).

V-D Dictionary Learning with BLOTLESS Update

This subsection compares dictionary learning performance for different dictionary update methods. The sparse coding algorithm is OMP. IterTLS in Section IV-A3 is used for BLOTLESS. Results for synthetic data are presented in Section V-D1 while Section V-D2 focuses on image denoising using real data.

V-D1 Synthetic Data

Fig. 7 and 8 compare the performance of dictionary learning using different dictionary update algorithms. Fig. 7 focus on the noise free case where $\bm{Y}=\bm{D}_{0}\bm{X}_{0}$ and Fig. 8 concerns with the noisy case where $\bm{Y}=\bm{D}_{0}\bm{X}_{0}+\bm{W}$ where $\bm{W}$ is the additive Gaussian noise matrix with i.i.d. entries and the signal-to-noise ratio (SNR) is set to 15dB. Both figures include the cases of complete and over-complete dictionaries. The results presented in Fig. 7 and 8 clearly indicate that dictionary learning based on BLOTLESS converges much faster and needs at least $1/3$ less training samples than other benchmark dictionary update methods.

In BLOTLESS update, blocks of the dictionary (sub-dictionaries) are updated sequentially. Fig. 9 compares the performance with different block sizes. Note that when each block contains only one dictionary item, the dictionary update problem is the same as that in K-SVD. Hence the performance of K-SVD is added in Fig. 9. Simulations suggest that the larger the dictionary blocks are, the faster the convergence is and the better performance is. The performance of BLOTLESS with block size one is slightly better than that of K-SVD. This is because the TLS step in IterTLS does not enforce the sparsity pattern and hence better accommodates errors in the estimated sparsity pattern.

V-D2 Real Data

We use the Olivetti Research Laboratory (ORL) face database [48] for dictionary learning and then use the learned dictionary for image denoising.

For dictionary learning, according to the simulation results in Section V-D1, $n=500$ samples of size $8\times 8$ patches from face images are enough for training a dictionary via BLOTLESS. The parameters used in dictionary learning are $m=64$, $l=128$, $\theta=10/128$, and $n_{\rm it}=20$. After learning a dictionary, image denoising [1] is performed using test images from the same dataset. The denoising results are shown in Table II, where four test images are used. In all these four tests, the BLOTLESS method outperforms all other algorithms, which is consistent with these simulations in V-D1.

VI Conclusion

This paper proposed a BLOTLESS algorithm for dictionary update. It divides the dictionary into sub-dictionaries, each of which is non-overcomplete. Then BLOTLESS updates a sub-dictionary and the corresponding sparse coefficients using least sqaures or total least squares approaches. Necessary conditions for unique recovery are identified and they hold with high probability when the number of training samples is larger than the derived bounds in Proposition 3. Simulations show that these bounds match the simulations well, and that BLOTLESS outperforms other benchmark algorithms. One future direction is to find sufficient bounds for unique recovery and their comparisons to the necessary bounds.

-A Proof of Proposition 3

The proof needs Hoeffding’s inequality[49] for Bernoulli random variables, stated below.

Lemma 1** (Hoeffding’s Inequality).**

For $N$ many identical Bernoulli random variables $\{X_{i},i=1,2,...N\}$. Each $X_{i}$ takes the value 1 with probability $p$ and 0 with probability $(1-p)$, then the following Hoeffding’s inequality holds

[TABLE]

where $\lambda>0$ is a constant number.

To derive $n_{1}$, we consider the case that the necessary condition 1 in Proposition 2 fails. That is,

[TABLE]

Based on Hoeffding’s inequality, the probability of this event is upper bounded by

[TABLE]

If this probability is smaller than $\epsilon$, it follows that

[TABLE]

The left hand side of the above inequality is quadratic in $n$, which after some elementary calculations leads to

[TABLE]

The derivation of $n_{2}$ is similar. Consider the probability that the necessary condition 2 in Proposition 2 fails:

[TABLE]

where the inequality in the third line follows from the union bound. After applying Hoeffding’s inequality and setting the upper bound less than $\epsilon$ we obtain

[TABLE]

It follows that

[TABLE]

To derive $n_{3}$, we define the following event

• •

$S_{i,i^{\prime}}$: For given $i\neq i^{\prime}$, $\exists j\in[n]$ such that $X_{i,j}=0$ and $X_{i^{\prime},j}\neq 0$.

Then the probability that the necessary condition 3 fails is given by

[TABLE]

where the inequality in the second line follows from the union bound. If we set this probability to be smaller than $\epsilon$, we obtain

[TABLE]