11institutetext: Algorithms and Complexity Group, TU Wien, Vienna, Austria
11email<EMAIL_ADDRESS>
# Layered Area-Proportional
Rectangle Contact Representations††thanks: We acknowledge funding by the
Austrian Science Fund (FWF) under grant P31119.

Martin Nöllenburg 11 Anaïs Villedieu 11 Jules Wulms 11

###### Abstract

We investigate two optimization problems on area-proportional rectangle
contact representations for layered, embedded planar graphs. The vertices are
represented as interior-disjoint unit-height rectangles of prescribed widths,
grouped in one row per layer, and each edge is ideally realized as a rectangle
contact of positive length. Such rectangle contact representations find
applications in semantic word or tag cloud visualizations, where a collection
of words is displayed such that pairs of semantically related words are close
to each other. In this paper, we want to maximize the number of realized
rectangle contacts or minimize the overall area of the rectangle contact
representation, while avoiding any false adjacencies. We present a network
flow model for area minimization, a linear-time algorithm for contact
maximization of two-layer graphs, and an ILP model for maximizing contacts of
$k$-layer graphs.

###### Keywords:

contact graphs layered planar graphs semantic word clouds

##  1 Introduction

Contact representations of planar graphs are a well-studied topic in graph
theory, graph drawing, and computational geometry [7, 10, 8]. Vertices are
represented by geometric objects, e.g., disks or polygons, and two objects
touch if and only if they are connected by an edge. They find many
applications, for instance in VLSI design [18], cartograms [13], or semantic
word clouds [17, 1].

Figure 1: Word cloud generated from the apnews.com frontpage by worditout.com
on the day of the certification of Joe Biden’s victory in the 2020 US
elections.

Word or tag clouds are popular visualizations that summarize textual
information in an aesthetically pleasing way. They show the main themes of a
text by displaying the most important keywords obtained from text analysis and
scale the word size to their frequency in the text. Word clouds became
widespread after the first automated generation tool “Wordle” was published in
2009 [14].

Word clouds with their different font sizes and words packed without semantic
context, such as the one shown in Fig. 1, have also received some criticism as
their audience sometimes fails at understanding the underlying data (while
enjoying their playful nature) [9]. For example, neighboring words that are
not semantically related can be misleading (see marked words in Fig. 1). As a
way to improve readability, _semantic_ word clouds have been introduced [6,
17, 1]. In semantic word clouds, an underlying edge-weighted graph indicates
the semantic relatedness of two words, whose positions are chosen such that
semantically related words are next to each other while unrelated words are
kept far apart.

Classic word clouds are often generated using forced-based approaches,
alongside with a spiral placement heuristic [14, 15, 16] that allows for a
very compact final layout. This method is powerful even when the rough
position of a word is dictated by an underlying map [4, 12]. Semantic word
clouds on the other hand have been approached with many different techniques,
e.g., force directed [6], seam-carving [17], and multidimensional scaling [2].
The problem has also been studied from a theoretical point of view, where an
edge of the semantic word graph is realized if the bounding boxes of two
related words properly touch; the realized edge weight is gained as profit.
Then the semantic word cloud problem can be phrased as the optimization
problem to maximize the total profit. Barth et al. [1] and later Bekos et al.
[3] gave several hardness and approximation results for this problem (and some
variations) on certain graph classes. The underlying geometric problem also
has links to more general contact graph representation problems, like
rectangular layouts [5] or cartograms [11].

In most of the literature about layered graphs, vertices are assigned to rows
without a predefined left-to-right order, yet this has interesting properties
in the context of word clouds. For instance, layered rectangle contact
representations are compact, assuming a good assignment they have an even
distribution of words and our eye naturally understands words grouped in rows
or tables. In this paper we study row-based contact graphs of unit-height but
arbitrary-width rectangles, which may represent the bounding boxes of words
with fixed font size.

Problem description. As input we take a layered graph $G=(V,E)$ on $L$ layers,
with an arbitrary number of vertices per layer. Each vertex $v_{i,j}\in V$ is
indexed by its layer $i\in[0,L-1]$ and its position $j$ within the layer:
$v_{i,j}$ is the $j^{th}$ vertex on the $i^{th}$ layer. The edge set $E$
consists of edges connecting each vertex $v_{i,j}$ to its neighbors
$v_{i,j-1}$ and $v_{i,j+1}$ on the same row (if they exist), and connections
between adjacent rows form an internally triangulated graph. We associate each
vertex with an axis-aligned unit-height rectangle $R_{i,j}$ with width
$w_{i,j}$, and $y$-coordinate $i$. We want to compute its $x$-position
$x_{i,j}$ given by the $x$-coordinate of its bottom left corner such that the
rectangles do not overlap except on their boundaries (see Fig. 2). Leaving
whitespace between two rectangles on the same layer is allowed and forms a
_gap_. Such a layout $\mathcal{R}$ is called a _representation_ of $G$. An
edge $(u,v)\in E$ is _realized_ in a representation $\mathcal{R}$ if
rectangles $R_{u}$ and $R_{v}$, representing vertices $u$ and $v$, intersect
along their boundaries for a positive length $\varepsilon>0$, which we denote
by $(R_{u},R_{v})\in\mathcal{R}$. If $R_{u}$ and $R_{v}$ are horizontally
adjacent we call the contact a _horizontal contact_.

Otherwise, the intersection is located along a horizontal boundary if $R_{u}$
and $R_{v}$ are on adjacent layers; these are called _vertical contacts_.
Contacts between rectangles whose vertices are not adjacent in $G$ are _false
adjacencies_. Such adjacencies can mislead a user to infer a link between
unrelated words, invalidating the representation. Within this model we study
two problem variations, _area minimization_ and _contact maximization_.

Figure 2: Partial drawing of a graph $G$, along with a representation
$\mathcal{R}$ of the visible vertices of $G$. Red edges are not realized, due
to the gray gap in $\mathcal{R}$.

For the area minimization problem the goal is to produce a representation
$\mathcal{R}$ that minimizes the total width of the gaps in $\mathcal{R}$. The
contact maximization problem asks to maximize the number of adjacencies
realized in $\mathcal{R}$, as specified by edge set $E$. For both optimization
criteria, false adjacencies are forbidden: otherwise a trivial gap-less
representation would always be a solution to the area minimization problem and
in the case of contact maximization, false adjacencies may reduce the number
of lost contacts with respect to a valid optimal solution as Fig. 3 shows. We
say a representation is _valid_ if it has no false adjacencies.

Figure 3: Allowing false adjacencies (red) could reduce lost adjacencies
(blue).

##  2 Area minimization

To solve the area minimization problem, we construct a flow network
$N=(G^{\prime}=(V^{\prime},E^{\prime});l;c;b;cost)$ for a given vertex-
weighted layered graph $G=(V,E)$, with edge capacity lower bound $l\colon
E^{\prime}\rightarrow\mathbb{R}^{+}_{0}$, edge capacity $c\colon
E^{\prime}\rightarrow\mathbb{R}^{+}_{0}$, vertex production/consumption
$b\colon V^{\prime}\rightarrow\mathbb{R}$ and cost function $cost\colon
E^{\prime}\rightarrow\mathbb{R}^{+}_{0}$. Each unit of cost will represent a
unit length gap and each unit of flow on an edge will represent a unit length
contact. To build the network we create two vertices $v^{a}$ and $v^{b}$ for
each rectangle, that respectively receive the flow from the lower layer and
output flow to the upper layer, and one for each potential gap, located
between each sequential pair of rectangles in the same layer. Every edge $e$
that ends on a gap vertex has $cost(e)=1$. We also add an edge $e$ between
$v^{a}$ and $v^{b}$ for each $R_{i,j}$ with $l(e)=c(e)=w_{i,j}$ and no cost to
ensure that rectangle nodes receive exactly as much flow as they are wide.

The intuition behind the network is that it represents a stack of layers
consisting of rectangles and gaps, with a maximum width of $w_{\max}\cdot K$,
$K$ being the maximum number of rectangles per layer, and $w_{\max}$ the width
of the widest rectangle. To facilitate this flow on all layers, there are
buffer vertices on both sides of each layer. Each rectangle is as wide as the
amount of flow its vertices $v^{a}$ and $v^{b}$ receive, and has contacts with
its upper and lower neighbors as wide as the flow on the edges representing
these contacts. Every vertex has edges to the layer above as far as its
rectangle is allowed to have contacts: a rectangle $R_{i,j}$ that has only one
upper neighbor, will have an edge to that neighbor, and to the gaps on that
neighbor’s right and left side. Any further edge would be to another rectangle
with which $R_{i,j}$ should not share a contact, and such edges would hence
result in false adjacencies. We picture the stack bottom-up, meaning that the
flow comes in at the bottom layer and exits from the top layer. A gap block
$g_{i,j}$ will reach as far left and right as its left and right neighbors in
the same row: if rectangle $R_{i,j}$ lies directly left of $g_{i,j}$, then the
furthest left upward neighbor of $R_{i,j}$ is the furthest left upward
neighbor of $g_{i,j}$. If $g_{i,j}$ could reach even further, then it would
essentially push $R_{i,j}$ into a false adjacency. The exact construction is
detailed in Appendix 0.A and sketched in Fig. 4.

Figure 4: Parts of a flow network: outgoing edges from the source (orange), a
rectangle (green) and a gap (blue); red buffer rectangles; gray edges have
cost 1.

###### Theorem 2.1

Given a graph $G=(V,E)$, the cost of a minimum-cost flow $f$ in $N$ equals the
minimum total gap length of any valid representation of $G$. An area-minimal
representation of $G$ is constructed from $f$ in polynomial time.

###### Proof

Given any graph $G=(V,E)$, the associated network $N$ has production equal to
its consumption $\sum_{v\in V^{\prime}}b(v)=b(s)+b(t)=0$. The source produces
$b(s)=w_{max}\cdot K$ flow, which is available to every vertex on layer $1$
(except $v^{b}$ vertices whose incoming neighbor is $v^{a}$ on the same
layer). Any vertex $v^{a}_{1,j}$ must receive $w_{1,j}$ units of flow, as its
only edge towards $v^{b}_{1,j}$ has capacity constraint $c=l=w_{1,j}$. Because
$w_{1,j}\leq w_{max}$ and since there are at most $K$ vertices $v^{a}_{1,j}$
on layer $1$, the capacities can be satisfied. Any edge in $E^{\prime}$ goes
from layer $i$ to $i+1$, except for $(v^{a}_{i,j},v^{b}_{i,j})$, but the exact
amount of flow that comes from layer $i-1$ into $v^{a}_{i,j}$ will go through
$v^{b}_{i,j}$ to layer $i+1$. Hence for the same reason as in layer $1$, there
is enough flow to satisfy the edge capacity constraints, while excess flow is
routed through (0-cost) buffers or (1-cost) gaps.

Since in this network only flow that goes into a gap vertex has a (non-zero)
cost, flow into gap vertices, and therefore also the total gap width, is
minimized. The minimum cost procedure finds this optimal flow $f$.

We construct a minimum-area representation $\mathcal{R}$ by placing rectangles
row by row: we leave buffers and gaps equal to the flow routed through the
corresponding vertices, and align all rows on the left. The total width of
each row, including buffers and gaps, is $w_{max}\cdot K$. Since the flow
through each rectangle is constant, the area of the buffers is maximized, to
minimize the area occupied by gaps. ∎

While this method minimizes the area occupied by the drawing it will not
always lead to the representation with the minimum bounding box. To minimize
the size of the bounding box, we propose to limit the amount of flow outgoing
from the source node and incoming into the sink node. If the chosen bound is
too small then the flow network will not be realizable. We can thus perform a
binary search between the width of the longest layer $W_{\max}$ as a lower
bound and $w_{\max}\cdot K$ as an upper bound. Assuming input widths have
integer values, this method would add a $\mathcal{O}(\log(w_{\max}\cdot
K-W_{\max}))$ factor to the flow runtime.

##  3 Maximization of realized contacts

In this section we propose algorithms that maximize the number of realized
contacts. We start with a linear-time algorithm for $L=2$, followed by an
integer linear programming model for $L>2$. The complexity for $L>2$ remains
open.

###  3.1 Linear-time algorithm for $L=2$

In this section we describe an algorithm $\mathcal{A}$ for the case where the
input has 2 layers. On a 2-layer graph, a vertex either has one neighbor in
the adjacent layer, or more than one. If a vertex has one neighbor we call it
a _T-vertex_ and if there are more neighbors, it is called a _fan_. A _block_
is a maximal sequence of consecutive rectangles in a layer $i$, for which each
horizontal contact is realized. A block from the $j$th until the $l$th vertex
of row $i$ is the sequence $(R_{i,j},\ldots,R_{i,l})$, where for each
$k\in[j,l-1]$ holds that $(R_{i,k},R_{i,k+1})\in\mathcal{R}$.

In a given 2-layer graph, there will always be a layer that starts on a fan,
while the opposite layer starts on a T-vertex. Assume without loss of
generality that $R_{0,0}$ is a fan, otherwise swap the two rows for the
duration of the algorithm. Algorithm $\mathcal{A}$ starts by placing
$R_{0,0}$, followed by all its neighbors on the adjacent layer, from left to
right, ending with $R_{1,j}$. Rectangle $R_{1,j}$ is again a fan, and the
process of placing all opposite-row neighbors, left to right, is repeated for
$R_{1,j}$ and every consecutive fan, as they are encountered. We call this
placement ordering $\prec$.

When we add a rectangle $R_{i}$ (fan or T-vertex), we always first attempt to
add it next to its horizontal predecessor, if possible (no false adjacency).
Though, if the horizontal predecessor is too far left, we place $R_{i}$ in the
leftmost allowed position. Let $R_{0}$ be the first rectangle in $\prec$,
which is placed on position $x_{0}$ by $\mathcal{A}$. Algorithm $\mathcal{A}$
then proceeds by adding $R_{1}$, representing a T-vertex in the opposite row.
Rectangle $R_{1}$, with width $w_{1}$, is placed leftmost, on coordinate
$x_{0}+\varepsilon-w_{1}$. We then proceed to add all rectangles corresponding
to other T-vertices of $R_{0}$ one by one, such that all horizontal contacts
are realized. Once a T-vertex $R_{i}$ cannot reach $R_{0}$, we store the
amount of contacts currently realized by $R_{0}$ as well as its position
$x_{0}$, and slide $R_{0}$ rightward, to the leftmost position
$x_{0}^{\prime}$ that allows a contact of $\varepsilon$ with $R_{i}$. Note
that, since we placed $R_{1}$ in the leftmost position that allowed a contact
of width $\varepsilon$ with $R_{0}$, we lose at least one contact by moving
$R_{0}$ rightward. If placing $R_{0}$ at $x_{0}^{\prime}$ ties the amount of
contacts of $x_{0}$, then we set $x_{0}\coloneqq x_{0}^{\prime}$. If
$x_{0}^{\prime}$ is strictly worse, then the representation is reset to having
$R_{0}$ at $x_{0}$. From that point on, every time we add a new rectangle, we
attempt this shift of the fan and update the position when we find a tie or
when we realize more contacts. We repeat this operation for each rectangle,
following the order $\prec$, always shifting the last encountered fan.

However, once we consider a fan $R_{f}$ that is not $R_{0}$, any sliding
operation will be attempted on the block containing $R_{f}$, rather than just
$R_{f}$. As before, we always shift the block to the leftmost position that
realizes the contact between $R_{f}$ and the newly placed rectangle. We
remember position $x_{f}$ that realizes most contacts, and favor the newest
position on a tie. In case moving the block containing $R_{f}$ leads to
strictly less contacts, we also try to move only $R_{f}$ instead. This starts
a new block containing just $R_{f}$. Below we sketch the proof for Theorem
3.1, the complete proof can be found in Appendix 0.B.

###### Theorem 3.1

Algorithm $\mathcal{A}$ computes a contact maximal valid representation with
contacts of length at least $\varepsilon$ for a given 2-layer graph $G$ in
linear time.

######  Proof (Sketch)

We show that during algorithm $\mathcal{A}$, the invariant holds that a
representation of the first $n$ rectangles in $\prec$ maximizes the number of
contacts.

We assume that the invariant holds after $\mathcal{A}$ placed $n-1$
rectangles, such that the current representation $\mathcal{R}^{*}$ is a
contact maximal representation of the first $n-1$ rectangles in $\prec$, and
realizes $k$ contacts. Algorithm $\mathcal{A}$ now adds the next rectangle
$R_{n}$. We show that the new representation is contact maximal.

Figure 5: Two configurations where $R_{l}$ is necessarily a fan (blue) or a
T-vertex (yellow). (a) If $R_{n}$ can achieve only a vertical contact, $R_{f}$
is a fan. (b) If $R_{n}$ can achieve only a horizontal contact and is a fan,
$R_{f}$ is a T-vertex.

We first prove that the maximum number of contacts that the new representation
can realize is $k+2$, since $R_{n}$ can achieve at most one vertical and one
horizontal contact. If these contacts happen naturally, when placing $R_{n}$
leftmost, then the invariant trivially holds. We therefore prove via a case
distinction that in all other cases $k+1$ adjacencies are optimal. We
distinguish between the contact that is achieved by $R_{n}$, either vertical
or horizontal. A sole vertical contact with fan (necessarily, forced by the
placement ordering) $R_{f}$ is achieved only if the horizontal predecessor
$R_{l}$ of $R_{n}$ is a fan, as shown in Fig. 5a. If $R_{n}$ is a T-vertex, a
single horizontal contact can arise only if $R_{f}$ is not moved to $R_{n}$,
as can be seen in Fig. 6. However, when $R_{n}$ is a fan, this requires
$R_{l}$ to be a T-vertex, see Fig. 5b. The invariant will therefore still be
true after $\mathcal{A}$ added all rectangles, producing a contact maximal
representation.

Algorithm $\mathcal{A}$ considers each rectangle either once, or its degree
many times, when a fan is shifted. As the input graph $G$ is planar,
$\mathcal{A}$ runs in linear time. ∎

Figure 6: Three configurations where T-vertex $R_{n}$ does not realize
vertical contacts with $R_{f}$ initially. We move $R_{f}$ and either (a) reset
if the number of contacts is strictly worse, or save when we find (b) a tie,
or (c) an increase in contacts.

###  3.2 ILP

To solve the contact maximization problem on $L>2$ layers we propose an ILP
formulation, which intuitively works as follows. We create a binary contact
variable $c(e)$ for each edge $e$ in the input graph. If a contact is not
realized, we set $c(e)=1$ to satisfy the position constraints, otherwise we
can set $c(e)=0$. To handle false adjacencies we add for each rectangle a
constraint on the first false contact that happens from the right and left on
the row above, if they exist. We use hard constraints on the rectangle
coordinates to prevent the false adjacencies. The objective is to minimize the
sum over all contact variables, under all these constraints, to maximize the
number of realized contacts in a solution. Additional details can be found in
Appendix 0.C.

$\text{minimize}\sum_{(v,v^{\prime})\in E}c(v,v^{\prime})$ (1)

The following inequalities ensure that there is no overlap between rectangles
on the same layer (2), and check whether the horizontal contact is realized
(3).

$\displaystyle x_{i,j}+w_{i,j}$ $\displaystyle\leq x_{i,j+1}$
$\displaystyle\forall(v_{i,j},v_{i,j+1})\in E$ (2) $\displaystyle x_{i,j+1}$
$\displaystyle\leq x_{i,j}+w_{i,j}+c(v_{i,j},v_{i,j+1})M$
$\displaystyle\forall(v_{i,j},v_{i,j+1})\in E$ (3)

The next inequalities verify that the contacts between rectangle $R_{i,j}$ and
all of its neighbors on layer $i+1$ are realized.

$\displaystyle x_{i+1,j^{\prime}}$ $\displaystyle\leq
x_{i,j}+w_{i,j}-\varepsilon+c(v_{i,j},v_{i+1,j^{\prime}})M$
$\displaystyle\forall e(v_{i,j},v_{i+1,j^{\prime}})\in E$ (4) $\displaystyle
x_{i,j}$ $\displaystyle\leq
x_{i+1,j^{\prime}}+w_{i+1,j^{\prime}}-\varepsilon+c(v_{i,j},v_{i+1,j^{\prime}})M$
$\displaystyle\forall e(v_{i,j},v_{i+1,j^{\prime}})\in E$ (5)

Finally, we model false adjacencies using pairs $(v_{i,j},v_{i+1,j^{\prime}})$
stored in sets $F_{L}$ (resp. $F_{R}$) that indicate the index of the first
block in row $i+1$ that is left (resp. right) of a neighbor of $R_{i,j}$, but
is not itself a neighbor of $R_{i,j}$.

$\displaystyle x_{i+1,j^{\prime}}+w_{i+1,j^{\prime}}$ $\displaystyle\leq
x_{i,j}$ $\displaystyle\forall(v_{i,j},v_{i+1,j^{\prime}})\in F_{L}$ (6)
$\displaystyle x_{i,j}+w_{i,j}$ $\displaystyle\leq x_{i+1,j^{\prime}}$
$\displaystyle\forall(v_{i,j},v_{i+1,j^{\prime}})\in F_{R}$ (7)

## References

  * [1] Barth, L., Fabrikant, S.I., Kobourov, S.G., Lubiw, A., Nöllenburg, M., Okamoto, Y., Pupyrev, S., Squarcella, C., Ueckerdt, T., Wolff, A.: Semantic word cloud representations: Hardness and approximation algorithms. In: Pardo, A., Viola, A. (eds.) Theoretical Informatics (LATIN’14). LNCS, vol. 8392, pp. 514–525. Springer (2014). https://doi.org/10.1007/978-3-642-54423-1_45
  * [2] Barth, L., Kobourov, S.G., Pupyrev, S.: Experimental comparison of semantic word clouds. In: Gudmundsson, J., Katajainen, J. (eds.) Experimental Algorithms (SEA’14). pp. 247–258. LNCS, Springer (2014). https://doi.org/10.1007/978-3-319-07959-2_21
  * [3] Bekos, M.A., van Dijk, T.C., Fink, M., Kindermann, P., Kobourov, S., Pupyrev, S., Spoerhase, J., Wolff, A.: Improved approximation algorithms for box contact representations. Algorithmica 77(3), 902–920 (2017). https://doi.org/10.1007/s00453-016-0121-3
  * [4] Buchin, K., Creemers, D., Lazzarotto, A., Speckmann, B., Wulms, J.: Geo word clouds. In: IEEE Pacific Visualization (PacificVis’16). pp. 144–151 (2016). https://doi.org/10.1109/PACIFICVIS.2016.7465262
  * [5] Buchsbaum, A.L., Gansner, E.R., Procopiuc, C.M., Venkatasubramanian, S.: Rectangular layouts and contact graphs. ACM Transactions on Algorithms 4(1), 8:1–8:28 (2008). https://doi.org/10.1145/1328911.1328919
  * [6] Cui, W., Wu, Y., Liu, S., Wei, F., Zhou, M.X., Qu, H.: Context preserving dynamic word cloud visualization. In: IEEE Pacific Visualization (PacificVis’10). pp. 121–128 (2010). https://doi.org/10.1109/PACIFICVIS.2010.5429600
  * [7] Felsner, S.: Rectangle and square representations of planar graphs. In: Pach, J. (ed.) Thirty Essays on Geometric Graph Theory, pp. 213–248. Springer (2013). https://doi.org/10.1007/978-1-4614-0110-0_12
  * [8] de Fraysseix, H., de Mendez, P.O., Rosenstiehl, P.: On triangle contact graphs. Combinatorics, Probability and Computing 3, 233–246 (1994). https://doi.org/10.1017/S0963548300001139
  * [9] Hearst, M.A., Pedersen, E., Patil, L., Lee, E., Laskowski, P., Franconeri, S.: An evaluation of semantically grouped word cloud designs. IEEE Transactions on Visualization and Computer Graphics 26(9), 2748–2761 (2020). https://doi.org/10.1109/TVCG.2019.2904683
  * [10] Koebe, P.: Kontaktprobleme der konformen Abbildung. Ber. Sächs. Akad. Wiss. Leipzig, Math.-Phys. Klasse 88, 141–164 (1936)
  * [11] van Kreveld, M., Speckmann, B.: On rectangular cartograms. Computational Geometry 37(3), 175–187 (2007). https://doi.org/10.1016/j.comgeo.2006.06.002
  * [12] Li, C., Dong, X., Yuan, X.: Metro-wordle: An interactive visualization for urban text distributions based on wordle. Visual Informatics 2(1), 50–59 (2018). https://doi.org/10.1016/j.visinf.2018.04.006
  * [13] Nusrat, S., Kobourov, S.: The state of the art in cartograms. Computer Graphics Forum 35(3), 619–642 (2016). https://doi.org/10.1111/cgf.12932
  * [14] Viegas, F.B., Wattenberg, M., Feinberg, J.: Participatory visualization with wordle. IEEE Transactions on Visualization and Computer Graphics 15(6), 1137–1144 (2009). https://doi.org/10.1109/TVCG.2009.171
  * [15] Wang, Y., Chu, X., Bao, C., Zhu, L., Deussen, O., Chen, B., Sedlmair, M.: Edwordle: Consistency-preserving word cloud editing. IEEE Transactions on Visualization and Computer Graphics 24(1), 647–656 (2018). https://doi.org/10.1109/TVCG.2017.2745859
  * [16] Wang, Y., Chu, X., Zhang, K., Bao, C., Li, X., Zhang, J., Fu, C.W., Hurter, C., Deussen, O., Lee, B.: Shapewordle: Tailoring wordles using shape-aware archimedean spirals. IEEE Transactions on Visualization and Computer Graphics 26(1), 991–1000 (2020). https://doi.org/10.1109/TVCG.2019.2934783
  * [17] Wu, Y., Provan, T., Wei, F., Liu, S., Ma, K.L.: Semantic-preserving word clouds by seam carving. Computer Graphics Forum 30(3), 741–750 (2011). https://doi.org/10.1111/j.1467-8659.2011.01923.x
  * [18] Yeap, K.H., Sarrafzadeh, M.: Floor-planning by graph dualization: 2-concave rectilinear modules. SIAM J. Comput. 22(3), 500–526 (1993). https://doi.org/10.1137/0222035

##  Appendix 0.A Flow construction

For each $v_{i,j}\in V$ we create two copies $v^{a}_{i,j}$ and $v^{b}_{i,j}$
in $V^{\prime}$; for each pair $v_{i,j},v_{i,j+1}$ we introduce a gap vertex
$g_{i,j}$, for each layer a left and right buffer vertex $l_{i}$ and $r_{i}$
and a global source and sink, $s$ and $t$. All the vertices $v$ have $b(v)=0$
except $b(s)=w_{max}\cdot K$ and $b(t)=-w_{max}\cdot K$ with $w_{max}$ the
width of the widest rectangle and $K$ the maximum number of rectangles per
layer.

Unless stated otherwise, each edge $e$ has $c(e)=\infty$, $l(e)=0$ and
$cost(e)=0$. The rectangle-rectangle edges for two vertically adjacent
vertices $u$ and $v$ are from the $u^{b}$ vertex on layer $i$ to the $v^{a}$
vertex on layer $i+1$. For each $v^{a}_{i,j},v^{b}_{i,j}$ pair in $G^{\prime}$
we add an edge $e_{i,j}$ from $v^{a}_{i,j}$ to $v^{b}_{i,j}$ with
$c(e_{i,j})=w_{i,j}$ and $l(e_{i,j})=w_{i,j}$. This ensures that a vertex must
receive exactly as much flow as the width of the rectangle it represents. We
also add edges going in and out of the buffer vertices: to minimize the gaps
it is preferable for the outer rectangles to get missing flow from the buffers
on the outside or to route excess flow to these buffers.

Edges from rectangle vertices to gap vertices are defined as follows. For a
vertex $v^{b}_{i,j}$, we add an edge to all vertices $v^{a}$ representing its
neighbors in $G$ on layer $i+1$, and we add edges to the gaps left and right
of those neighbors. All edges to gaps have $cost=1$, to penalize the creation
of gaps. Lastly we add edges from gaps to other vertices. For a gap vertex
$g_{i,j}$ we look at its left and right neighbors $v^{b}_{i,j}$ and
$v^{b}_{i,j+1}$ respectively. We add edges from $g_{i,j}$ to all rectangles
and gaps that have an incoming edge from $v^{b}_{i,j}$ and $v^{b}_{i,j+1}$,
and assign $cost=0$ for edges into rectangles, and $cost=1$ for edges into
gaps. We do not penalize the outgoing flow of a gap, since we already count
the incoming flow.

From this network we can easily deduce that the minimum cost flow corresponds
to the solution that minimizes total gap length: only gap vertices have a
cost, hence the flow avoids these vertices whenever possible. We can construct
the network and then solve the minimum flow problem in polynomial time. In the
obtained solution, the flow values found on each edge should represent the
length of the overlap between these elements, which allow us to construct the
corresponding representation. However, in practice the relationship between
units of flow and contact length is not always direct. Indeed, consider a 4
vertex configuration where $a$ lies bottom left, $b$ bottom right, $c$ top
left and $d$ top right, and both $a$ and $b$ have an edge to both $c$ and $d$.
Note that in such a case $a$ and $d$, or $c$ and $b$, would be gaps. It is
possible that $a$ favors sending its flow to $d$ and $b$ to $c$, which is
impossible to represent as contacts in a configuration of rectangles.

Once the flow has been computed, we can locally swap the required flow between
the crossing edges to resolve those _crossing patterns_.

##  Appendix 0.B Linear-time algorithm – Omitted proof

###### Theorem 3.1

Algorithm $\mathcal{A}$ computes a contact maximal valid representation with
contacts of length at least $\varepsilon$ for a given 2-layer graph $G$ in
linear time.

###### Proof

We show that during the placement of rectangles by algorithm $\mathcal{A}$,
the invariant holds that a representation computed for the first $n$
rectangles in $\prec$ achieves a maximum number of contacts. First observe
that, we start by placing $R_{0}$ and $R_{1}$, such that their only contact is
realized. The invariant therefore holds at the start.

Now assume that the invariant holds after $\mathcal{A}$ placed $n-1$
rectangles, such that the current representation $\mathcal{R}^{*}$ is a
contact maximal representation of the first $n-1$ rectangles in $\prec$, and
realizes $k$ contacts. Algorithm $\mathcal{A}$ now adds the next rectangle
$R_{n}$. We show that the new representation is contact maximal.

The maximum number of contacts that the new representation can realize is
$k+2$, because rectangle $R_{n}$ can achieve at most one vertical and one
horizontal contact. If there was a representation of the $n$ rectangles that
realized $k^{\prime}>k+2$ contacts, then removing $R_{n}$ would leave
$k^{\prime}-2>k$ contacts, which contradicts our assumption that
$\mathcal{R}^{*}$ is optimal. If $R_{n}$ naturally realizes its vertical and
horizontal contacts, the representation realizes $k+2$ contacts. There is no
configuration where we cannot realize any contact when adding rectangle
$R_{n}$ to $\mathcal{R}^{*}$, since one of the rows extends further, and hence
$R_{n}$ either can achieve its horizontal contact, or the opposing row extends
at least $\varepsilon$ and the vertical contact can be realized. We therefore
consider only the cases where adding $R_{n}$ results in a representation with
$k+1$ contacts. Throughout the rest of the proof, we refer to the vertical fan
neighbor of $R_{n}$, that has already been placed, as $R_{f}$. $R_{f}$
necessarily lies in the opposite layer. The only rectangles that will be moved
when placing $R_{n}$ are the rectangles in the block containing $R_{f}$. Any
other rectangles and adjacencies are untouched. So it is sufficient to count
the adjcencies gained and lost by the block containing $R_{f}$ and $R_{n}$ to
confirm that the representation of $n$ rectangles is optimal.

  * •

If $R_{n}$ realizes only its vertical contact with $R_{f}$, then rectangle
$R_{l}$, neighboring $R_{n}$ on the left, is necessarily a fan (see Fig. 7a)).
Assume for contradiction that $R_{l}$ is a T-vertex. This implies that $R_{f}$
is $R_{l}$’s only allowed vertical contact. Since $R_{n}$ achieves contact
with $R_{f}$ and $R_{f}$ is leftmost, $R_{l}$ must also be in contact with
$R_{f}$. As a result the horizontal contact with $R_{l}$ always happens,
contradicting our assumption that $R_{n}$ has only a vertical contact. As a
consequence, $R_{l}$ is a fan, and the position for a fan is the position that
maximizes contacts, favoring newer positions (or positions more to the right)
for ties. An alternative to $\mathcal{R}^{*}$ with $R_{l}$ in contact with
$R_{f}$ must therefore achieve less than $k$ contacts. Additionally, $R_{n}$
cannot be slid left of $R_{f}$ even though there is a gap towards $R_{l}$, as
the rectangles left of $R_{f}$ are not part in the neighborhood of $R_{n}$.
Thus, sliding $R_{n}$ further left would cause false adjacencies. There is no
contact maximal representation that preserves the contact between $R_{l}$ and
$R_{n}$, and hence $k+1$ contacts is maximal.

  * •

If $R_{n}$ realizes only its horizontal contact, then it fails the contact
with $R_{f}$. We distinguish between $R_{n}$ being a fan or a T-vertex.

    * –

If $R_{n}$ is a T-vertex, then $R_{f}$ is too far to the left, and the
algorithm will try to create representations where $R_{f}$ is in the leftmost
position that realizes a contact with $R_{n}$: either by moving the block
containing $R_{f}$, or alternatively moving just $R_{f}$.

Figure 7: Two configurations where $R_{l}$ is necessarily a fan (blue) or a
T-vertex (yellow). (a) If $R_{n}$ can only achieve a vertical contact, $R_{f}$
is a fan. (b) If $R_{n}$ can only achieve a horizontal contact and is a fan,
$R_{f}$ is a T-vertex.

      * *

If $R_{f}$ remains in its original position, which does not achieve the
vertical contact with $R_{n}$, it means that the leftmost position of $R_{f}$,
that achieves the contact with $R_{n}$ realizes strictly less contacts than
the original position of $R_{f}$ in $\mathcal{R}^{*}$. As a result, there is
no contact maximal representation that preserves the contact between $R_{f}$
and $R_{n}$, and thus $k+1$ contacts is optimal (see Fig. 8a).

      * *

When the block containing $R_{f}$ is shifted to the right to create the
alternative representation that ties the number of contacts of
$\mathcal{R}^{*}$, at most one vertical contact must have been lost, to gain
the vertical one. Otherwise, $R_{f}$ can give up its horizontal contact to
achieve the same result (see Fig 8b). Thus, while $R_{n}$ realises 2 contacts,
$R_{f}$ lost one and the representation realizes $k+1$ contact. This
alternative representation is optimal and used by $\mathcal{A}$ instead of the
representation where $R_{f}$ is not moved rightward.

      * *

If the block containing $R_{f}$ does not lose contacts when shifted, then the
overall representation will gain the contact between $R_{n}$ and $R_{f}$
contact. This configuration is optimal because it achieves $k+2$ contacts (see
Fig. 8c).

    * –

If $R_{n}$ is a fan, then rectangle $R_{l}$, again neighboring $R_{n}$ on the
left, must be a T-vertex, since $R_{f}$ cannot be an empty fan (see Fig. 7b).
Assume for contradiction that $R_{l}$ is a fan, then $\mathcal{A}$ would try
to realize the vertical contacts between $R_{f}$ and the fans $R_{l}$ and
$R_{n}$ on the other row. In this case $R_{f}$ would necessarily be an empty
fan, and hence $\mathcal{A}$ will try to sacrifice the contact with the left
neighbor of $R_{f}$, to gain the vertical contact between $R_{f}$ and $R_{n}$.
As there are no T-vertices for $R_{f}$ to lose contacts with, $\mathcal{A}$
will always find at least a tie between the configuration that has the
vertical contact between $R_{f}$ and $R_{n}$, and $\mathcal{R}^{*}$, and hence
prefer the new configuration. This contradicts our assumption that $R_{n}$
realizes only its horizontal contact, and hence $R_{l}$ must be a T-vertex.
When $R_{n}$ is added, an alternative representation is created with $R_{f}$
shifted to the leftmost position that realizes the contact with $R_{n}$. Since
this configuration is not chosen by $\mathcal{A}$, the configuration realizes
less than $k+1$ contacts, and therefore $k+1$ contacts is optimal.

Thus, after adding $R_{n}$ to $\mathcal{R}^{*}$, algorithm $\mathcal{A}$
produces a representation that realizes $k+2$ contacts, if it is possible, and
$k+1$ otherwise, which is optimal. The invariant will therefore still be true
after $\mathcal{A}$ added all rectangles, producing a contact maximal
representation of our input graph.

The algorithm considers each rectangle either only once, or as many times as
its degree when an alternative shifted representation is created. As the input
graph is planar, the algorithm runs in linear time. ∎

Figure 8: Three configurations where T-vertex $R_{n}$ does not realize
vertical contacts with $R_{f}$ initially. We move $R_{f}$ and either (a) reset
if the number of contacts is strictly worse, or save when we find (b) a tie,
or (c) an increase in contacts.

As one may already have realized, in some cases, when looking at two
sequential fans that are not in contact, two T-vertices will be in a point
contact. This happens when the T-vertices cannot overlap without creating a
false adjacency, as we show in see Fig 9.

Figure 9: Edge case where the prefered position of the fan vertices (red and
blue) causes a point contact on fan vertices (orange and green)

##  Appendix 0.C ILP model

In this section we present an integer linear programming formulation $M$,
which can be used to find optimal solutions to the contact maximization
problem. In this formulation, the constants are the width $w_{i,j}$ of each
rectangle, the list $E$ of contacts, the minimal length $\varepsilon$ of a
contact and a large integer $M$ for indicator type constraints. Additionally,
we have the lists $F_{L}$ and $F_{R}$ of false adjacencies, defined as
follows. Consider the rectangle representing vertex $v_{i,j}$ and its
neighbors $N=\\{v_{i+1,k},\ldots,v_{i+1,l}\\}$ on the row above. For each such
vertex $v_{i,j}$, we store the pairs $(v_{i,j},v_{i+1,k-1})\in F_{L}$ and
$(v_{i,j},v_{i+1,l+1})\in F_{R}$, if $v_{i+1,k-1}$ and $v_{i+1,l+1}$ exist.
These pairs represent the false adjacencies of $R_{i,j}$ with the rightmost
rectangle left of $N$ and the leftmost rectangle right of $N$, respectively.
Note that preventing these two false adjacencies will prevent all false
adjacencies of $R_{i,j}$ with rectangles in row $i+1$, as long as the order of
the rectangles in both rows is correctly maintained.

The variables of the ILP are the following. We use $x_{i,j}$ as variables for
the position of each rectangle, representing the bottom left corner of the
rectangle, and Boolean variables $c(v,v^{\prime})$ which indicate that a
contact is not realized between $v$ and $v^{\prime}$. We start by stating the
optimization function, which minimizes the sum of all $c$ variables. In
practice, this means we maximize realized adjacencies.

$\text{minimize}\sum_{(v,v^{\prime})\in E}c(v,v^{\prime})$ (8)

The following inequalities ensure firstly, that there is no overlap between
rectangles on the same layer, and secondly, check whether the horizontal
contact is realized: if $x_{i,j}$ is too small, then $c$ must be set to $1$,
and hence the represented rectangle is too far left to have the contact
naturally.

$\displaystyle x_{i,j}+w_{i,j}$ $\displaystyle\leq x_{i,j+1}$
$\displaystyle\forall(v_{i,j},v_{i,j+1})\in E$ (9) $\displaystyle x_{i,j+1}$
$\displaystyle\leq x_{i,j}+w_{i,j}+c(v_{i,j},v_{i,j+1})M$
$\displaystyle\forall(v_{i,j},v_{i,j+1})\in E$ (10)

The next inequalities ensure that the contacts between rectangle $R_{i,j}$ and
all of its neighbors on layer $i+1$ are realized, such that the left side of a
neighbor is left of the right side of $R_{i,j}$, or symmetrically, the right
side of a neighbor is right of the left side of $R_{i,j}$. Again, if we cannot
satisfy the inequality with $c=0$, then it is set to 1, and since $M$ is a
large value, the inequalities are trivially satisfied.

$\displaystyle x_{i+1,j^{\prime}}$ $\displaystyle\leq
x_{i,j}+w_{i,j}-\varepsilon+c(v_{i,j},v_{i+1,j^{\prime}})M$
$\displaystyle\forall e(v_{i,j},v_{i+1,j^{\prime}})\in E$ (11) $\displaystyle
x_{i,j}$ $\displaystyle\leq
x_{i+1,j^{\prime}}+w_{i+1,j^{\prime}}-\varepsilon+c(v_{i,j},v_{i+1,j^{\prime}})M$
$\displaystyle\forall e(v_{i,j},v_{i+1,j^{\prime}})\in E$ (12)

Finally, we model false adjacencies using the pairs stored in $F_{L}$ and
$F_{R}$. For a pair in $F_{L}$, the rectangle $R_{i+1,j^{\prime}}$, with which
$R_{i,j}$ has a false adjacency, should stay left of $R_{i,j}$. This forces
the left side of $R_{i,j}$ to be right of the right side of
$R_{i+1,j^{\prime}}$. Symmetrically, pairs in $F_{R}$ prevent false
adjacencies when $R_{i+1,j^{\prime}}$ is to the right of $R_{i,j}$.

$\displaystyle x_{i+1,j^{\prime}}+w_{i+1,j^{\prime}}$ $\displaystyle\leq
x_{i,j}$ $\displaystyle\forall(v_{i,j},v_{i+1,j^{\prime}})\in F_{L}$ (13)
$\displaystyle x_{i,j}+w_{i,j}$ $\displaystyle\leq x_{i+1,j^{\prime}}$
$\displaystyle\forall(v_{i,j},v_{i+1,j^{\prime}})\in F_{R}$ (14)

###### Theorem 0.C.1

Solving ILP model $M$ optimally, results in an optimal solution for the
contact maximization problem.

###### Proof

First, we show that whenever a variable $c$ is set to zero, then a contact is
realized. Let us assume that $c=0$, then one of the following two cases
applies.

  * •

If $c(v_{i,j},v_{i,j+1})=0$, then the horizontal contact between $R_{i,j}$ and
$R_{i,j+1}$ is realized. Since $x_{i+1,j}\leq x_{i,j}+w_{i,j}$ and
$x_{i,j}+w_{i,j}\leq x_{i,j+1}$, then $x_{i,j}+w_{i,j}=x_{i,j+1}$. Thus the
coordinate of the right side of $R_{i,j}$ lies on the left side of
$R_{i,j+1}$, and there is a horizontal contact.

  * •

If $c(v_{i,j},v_{i+1,j^{\prime}})=0$, then the vertical contact between
$R_{i,j}$ and $R_{i+1,j^{\prime}}$ is realized. We have both
$x_{i+1,j^{\prime}}\leq x_{i,j}+w_{i,j}-\varepsilon$ which means that the left
side of $R_{i,j}$ lies left of the right side of $R_{i+1,j^{\prime}}$ and
$x_{i,j}\leq x_{i+1,j^{\prime}}+w_{i+1,j^{\prime}}$ which means that the right
side of $R_{i,j}$ lies right of the left side of $R_{i+1,j^{\prime}}$.

Finally, since we minimize the sum of the $c$ values, as many as possible are
set to $0$. Each $c$ value set to zero corresponds to a realized adjacency,
and thus the number or contacts in the resulting representation is maximized.
∎


# Reduced Kernel Dictionary Learning

Denis C. ILIE-ABLACHIM
Faculty of Automatic Control and Computers
University Politehnica of Bucharest
<EMAIL_ADDRESS>Bogdan DUMITRESCU
This work was supported by grants of the Romanian Ministry of Research,
Innovation and Digitization, CCCDI - UEFISCDI, project numbers PN-
III-P2-2.1-PED-2019-3248 and PN-III-P4-PCE-2021-0154, within PNCDI III.
Faculty of Automatic Control and Computers
University Politehnica of Bucharest
<EMAIL_ADDRESS>
###### Abstract

In this paper we present new algorithms for training reduced-size nonlinear
representations in the Kernel Dictionary Learning (KDL) problem. Standard KDL
has the drawback of a large size of the kernel matrix when the data set is
large. There are several ways of reducing the kernel size, notably Nyström
sampling. We propose here a method more in the spirit of dictionary learning,
where the kernel vectors are obtained with a trained sparse representation of
the input signals. Moreover, we optimize directly the kernel vectors in the
KDL process, using gradient descent steps. We show with three data sets that
our algorithms are able to provide better representations, despite using a
small number of kernel vectors, and also decrease the execution time with
respect to KDL.

##  I Introduction

Dictionary Learning (DL) is a representation learning method that aims to find
a sparse representation for a set of signals, $\bm{Y}$, represented as a
matrix with $N$ columns (signals) of size $m$. The representation is achieved
by computing a dictionary $\bm{D}$ of size $m\times n$ and a sparse
representation $\bm{X}$ of size $n\times N$ such that a good approximation
$\bm{Y}\approx\bm{DX}$ is obtained. Most applications with dictionary learning
are in image denoising, inpainting, signal reconstruction, clustering or
classification.

In this paper we present new methods of dictionary learning that produce
sparse representations in both linear and nonlinear spaces, starting from the
Kernel Dictionary Learning (KDL) idea [1, 2]. In a first approach, this is
done in two stages; a linear representation is built and the resulting
optimized dictionary is used unchanged in the nonlinear space as kernel
vectors during the training procedure. In the second approach, the kernel
vectors are optimized alongside with the nonlinear representation. The main
advantage of this method is the use of a reduced matrix, $\bm{D}$, containing
the kernel vectors, which is also the dictionary in a standard DL problem. By
the use of the reduced matrix, the built-in kernel is smaller and thus the
problem complexity is reduced.

As KDL has a high complexity when $N$ is large (and so the kernel matrix is
large), solutions have been adopted from other problems where kernels appear.
In Large-Scale Kernel Machines, various kernel enhancement or resizing
strategies have been used, such as Nyström Sampling [3], [4] or Random Fourier
Features (RFF) [5]. The first one computes a rank $\hat{m}$ approximation
$\hat{\bm{K}}$ of the kernel matrix $\bm{K}$. The whole procedure consists in
approximating the nonlinear mapping function $\varphi(\bm{Y})$ with a matrix
$\hat{\bm{Y}}$, containing a compressed version of the original signals.
Compared to the original problem, the new problem no longer requires a high
computational cost. The Random Fourier Features method proposes to map the
input signals to a randomized low-dimensional feature space. More exactly, the
inner product, used in the kernel trick, is replaced with a randomized map
$k(\bm{x},\bm{y})=\varphi(\bm{x})^{\top}\varphi(\bm{y})=z(\bm{x})^{\top}z(\bm{y})$,
where $z:\mathbb{R}^{m}\mapsto\mathbb{R}^{\hat{m}}$ and $\hat{m}\ll m$. Both
methods enable the use of fast linear methods, which will further use the
resulting features. The use of reduced kernels dates from Support Vector
Machine times, significant contributions being [6, 7], with applications in
classification [8]. Other learning methods where a reduced kernel appears can
be found in [9]. In these early works, kernel vectors are usually selected
(randomly or with some heuristic) as a subset of input signals. Finally, there
are KDL substitutes like [10], where the nonlinear transformation is performed
by a neural encoder-decoder, with standard DL on the encoded signals.

The paper is organized as follows. In Section II we review the standard DL and
KDL problems and the most common algorithm for solving them. In Section III we
present our own contribution, named the Reduced Kernel Dictionary Learning
(RKDL) problem, under three different scenarios, in which the kernel vectors
are: a) the result of a DL problem, solved a priori; b) optimized together
with the other KDL variables, using gradient descent, with an objective that
is directly related to that of the KDL problem; c) optimized like before, but
using a mixed objective that combines the nonlinear representation of the
signals with the linear representation of the kernel vectors. The algorithms
for b) and c) are new and their representation error is smaller. Section IV
contains the experimental results, obtained on three public data sets: Digits,
MNIST [11] and CIFAR-10 [12], under the three proposed scenarios.

##  II Kernel Dictionary Learning

The DL problem is formulated as following

$\begin{array}[]{ll}\displaystyle\min_{\bm{D},\bm{X}}&\|\bm{Y}-\bm{D}\bm{X}\|_{F}^{2}\\\
\text{ s.t. }&\left\|\bm{x}_{\ell}\right\|_{0}\leq s_{x},\ell=1:N\\\
&\left\|\bm{d}_{j}\right\|=1,j=1:n,\end{array}$ (1)

where $\left\|\cdot\right\|_{0}$ represents the $0$-pseudo-norm, $s_{x}$ is
the sparsity level and $\bm{d}_{j}$ is a column (named also atom) of $\bm{D}$.

The standard dictionary learning problem can be solved by using simple
strategies. In order to overcome the nonconvexity and the huge dimension of
the problem, the most usual optimization procedure iterates two basic steps
and is also known as DL by Alternate Optimization. In this way, the problem is
divided in two subproblems: sparse coding and dictionary update. By
alternating these two stages for a given number of iterations, good local
solutions can be obtained. A simple iteration consists of computing the sparse
representation $\bm{X}$, using Orthogonal Matching Pursuit (OMP) [13] while
the dictionary $\bm{D}$ is fixed, and then updating the dictionary atoms
successively while the sparse representation is fixed. There are several
methods of dictionary learning [14], but the one of interest to us is AK-SVD
[15], [16] due to its low complexity and good performance.

In the DL problem, the input data are modeled through a linear representation,
which in some cases may be seen as a limitation. In order to overcome this
drawback, kernel representations can be used for a better quantification of
similarities or differences between input vectors.

The kernel representation is an extension to nonlinearity. We do this by
associating to a data vector $\bm{y}\in\mathbb{R}^{m}$ a feature vector
$\varphi(\bm{y})\in\mathbb{R}^{\tilde{m}}$, where
$\varphi:\mathbb{R}^{m}\rightarrow\mathbb{R}^{\tilde{m}}$ is a nonlinear
function. Typically, Mercer kernels are used, which can be expressed as a
scalar product of feature vector functions
$k(\bm{x},\bm{y})={\varphi}(\bm{x})^{\top}{\varphi}(\bm{y})$. In the last
form, the scalar product can be replaced with a specific function definition,
such as radial basis function (RBF)
$k(\bm{x},\bm{y})=\exp(-\frac{\|\bm{x}-\bm{y}\|^{2}_{2}}{2\sigma^{2}})$ or
polynomial kernel $k(\bm{x},\bm{y})=(\bm{x}^{\top}\bm{y}+\alpha)^{\beta}$.

The Kernel Dictionary Learning (KDL) [1], [2] problem is

$\begin{array}[]{ll}\displaystyle\min_{\bm{A},\bm{Z}}&\|\varphi(\bm{Y})-\varphi(\bm{Y})\bm{A}\bm{Z}\|_{F}^{2}\\\
\text{ s.t. }&\left\|\bm{z}_{\ell}\right\|_{0}\leq s_{z},\ell=1:N\\\
&\left\|\varphi(\bm{Y})\bm{a}_{j}\right\|=1,j=1:n,\end{array}$ (2)

where $\varphi(\bm{Y})$ represents the nonlinear extension of data and
$\varphi(\bm{Y})\bm{A}$ is the kernel dictionary, where $\bm{A}$ is the
coefficients matrix of size $N\times p$. Depending on the used data set, the
problem can be difficult to solve due the large kernel matrix
$\bm{K}_{YY}=\varphi(\bm{Y})^{\top}\varphi(\bm{Y})$ that results from the
trace form of the objective function. In this case, the problems with large
data sets can involve large volume of memory and long execution times. The KDL
problem can also be solved by alternate optimization. The sparse
representation is computed according to the Kernel OMP algorithm [1]. The
columns of the matrix $\bm{A}$ are sequentially updated with an algorithm
inspired by AK-SVD, while the representation matrix $\bm{X}$ is fixed.

Both algorithms, AK-SVD and Kernel AK-SVD, alongside with the sparse
representation computation are presented in [14].

##  III Reduced Kernel Dictionary Learning

Nonlinear space can extend the horizon of data representation. However, KDL
has disadvantages when the number of available signals is large. The size of
the kernel increases in proportion to the size of the data. Thus, the problem
becomes more complex from a numerical point of view.

In order to overcome this limitations we use a reduced space on which the
kernel matrix is built. This strategy is implemented by using as kernel
vectors not the full set of signals, $\bm{Y}$, but a smaller set of vectors,
$\bm{D}$, trained with DL as a dictionary for linear representations, thus
replacing (2) with

$\begin{array}[]{ll}\displaystyle\min_{\bm{A},\bm{Z}}&\|\varphi(\bm{Y})-\varphi(\bm{D})\bm{A}\bm{Z}\|_{F}^{2}\\\
\text{ s.t. }&\left\|\bm{z}_{\ell}\right\|_{0}\leq s_{z},\ell=1:N\\\
&\left\|\varphi(\bm{D})\bm{a}_{j}\right\|=1,j=1:n.\end{array}$ (3)

This problem has two advantages that can be used as needed. If the number of
training signals is very large, a dictionary with a much smaller number of
atoms can be used. On the other hand, we have problems where the number of
training signals is small or the given signals are not representative enough
for the representation problem. In this case it is recommended to use large
dictionaries. In our work the case of interest is the first one.

###  III-A Standard Reduced Kernel Dictionary Learning

A first approach to (3) consists of solving it in two steps. In the first step
we design $\bm{D}$ by solving a DL problem and thus obtaining a trained
dictionary. We call this method RKDL-D, due the use of matrix $\bm{D}$; it was
introduced in [17]. Here is a brief reminder of the atom update step of
RKDL-D. (Sparse representation can be easily derived from Kernel OMP.)

Expressing the objective of (3) in its trace form and isolating the current
atom $a_{j}$, we can write

$\begin{split}\operatorname{Tr}\left[\left(\varphi^{\top}(\bm{Y})-\sum_{i\neq
j}\bm{z}_{i}\bm{a}_{i}^{\top}\varphi^{\top}(\bm{D})-\bm{z}_{j}\bm{a}_{j}^{\top}\varphi^{\top}(\bm{D})\right)\right.\\\
\left.\left(\varphi(\bm{Y})-\varphi(\bm{D})\sum_{i\neq
j}\bm{a}_{i}\bm{z}_{i}^{\top}-\varphi(\bm{D})\bm{a}_{j}\bm{z}_{j}^{\top}\right)\right].\end{split}$

We compute the partial derivatives with respect to the current atom
$\bm{a}_{j}$

$\frac{\partial\operatorname{Tr}(\cdot)}{\partial\bm{a}_{j}}=2\|\bm{z}_{j}\|^{2}\bm{K}_{DD}\bm{a}_{j}+2\bm{K}_{DD}\bm{R}\bm{z}_{j}-2\bm{K}_{YD}\bm{z}_{j}\\\
$

and the current sparse representation vector $\bm{z_{j}}$

$\frac{\partial\operatorname{Tr}(\cdot)}{\partial\bm{z}_{j}}=2\bm{a}_{j}^{\top}\bm{K}_{DD}\bm{a}_{j}\bm{z}_{j}+2\bm{R}^{\top}\bm{K}_{DD}\bm{a}_{j}-2\bm{K}_{YD}\bm{a}_{j}.$

where we have used the notations
$\bm{K}_{YD}=K(\bm{Y},\bm{D})=\varphi(\bm{Y})^{\top}\varphi(\bm{D})$,
$\bm{K}_{DD}=K(\bm{D},\bm{D})$ and $\displaystyle\bm{R}=\sum_{i\neq
j}\bm{a}_{i}\bm{z}_{i}^{\top}$.

Setting these derivatives to zero, we obtain the optimal atom (for fixed
representation)

$\bm{a}_{j}=\left(\|\bm{z}_{j}\|_{2}^{2}\bm{K}_{DD}\right)^{-1}(\bm{K}_{YD}+\bm{K}_{DD}\bm{R})\bm{z}_{j}$
(4)

and optimal representation (for fixed atom)

$\bm{z}_{j}=\left(\bm{K}_{YD}-\bm{R}^{\top}\bm{K}_{DD}\right)\bm{a}_{j}.$ (5)

The RKDL-D procedure of updating atoms is summarized in Algorithm 1. For
brevity, we did not use special notations, but only the signals from $\bm{Y}$
where $\bm{a}_{j}$ appears in the representation are involved in the
computation.

Data: complementary kernel matrix $\bm{K}_{DD}\in\mathbb{R}^{p\times p}$

partial kernel matrix $\bm{K}_{YD}\in\mathbb{R}^{N\times p}$

current kernel dictionary $\bm{A}\in\mathbb{R}^{N\times n}$

representation matrix $\bm{Z}\in\mathbb{R}^{n\times N}$

Result: updated kernel dictionary $\bm{A}$, representation $\bm{Z}$

1 Compute sum
$\displaystyle\bm{S}=\sum_{i=1}^{n}\bm{Z}_{i}^{\top}\bm{a}_{i}^{\top}$

2 for  _$j=1$ to $n$_ do

3 Modify sum: $\bm{R}=\bm{S}-\bm{Z}_{j}\bm{a}_{j}^{\top}$

4 Update atom:
$\bm{a}_{j}=\left(\|\bm{z}\|_{2}^{2}\bm{K}_{DD}\right)^{-1}(\bm{K}_{YD}^{\top}+\bm{K}_{DD}\bm{R})\bm{Z}_{j}$

5 Normalize atom:
$\bm{a}_{j}\leftarrow\bm{a}_{j}/\left(\bm{a}_{j}^{\top}\bm{K}_{DD}\bm{a}_{j}\right)^{\frac{1}{2}}$

6 Update representation:
$\bm{Z}_{j}^{\top}\leftarrow(\bm{K}_{YD}-\bm{R}\bm{K}_{DD})\bm{a}_{j}$

7 Recompute error: $\bm{S}=\bm{R}+\bm{Z}_{j}\bm{a}_{j}^{\top}$

8

Algorithm 1 RKDL-D – update step

###  III-B Optimized Reduced Kernel Dictionary Learning

The improvement that we propose here is to update the dictionary $\bm{D}$,
containing the kernel vectors, during the nonlinear optimization procedure. We
keep the idea of alternate optimization. The matrices $\bm{Z}$, $\bm{A}$ and
$\bm{D}$ are updated successively. As above, for $\bm{Z}$ we use Kernel OMP
and the atoms of $\bm{A}$ are updated as described by Algorithm 1. Updating
the dictionary $\bm{D}$ must be done with a different procedure, detailed
below. Since the dictionary $\bm{D}$ is updated together with the nonlinear
representation, we call the resulting method Optimized Reduced Kernel
Dictionary Learning (ORKDL-D).

In order to solve the optimization problem for $\bm{D}$, we update each column
$\bm{d_{j}}$ independently, by the use of the trace form of the objective
function of (3):

$\operatorname{Tr}\left[\bm{K}_{YY}-2\bm{K}_{YD}\bm{A}\bm{Z}+\bm{Z}^{\top}\bm{A}^{\top}\bm{K}_{DD}\bm{A}\bm{Z}\right].$

We compute the partial derivatives with respect to the $i$th element of the
current column $\bm{d}_{j}$, for both nonlinear terms

$\begin{array}[]{rcl}\frac{\partial\operatorname{Tr}\left[\bm{K}_{YD}\bm{A}\bm{Z}\right]}{\partial\bm{d_{ij}}}&=&\operatorname{Tr}\left[\left(\frac{\partial\operatorname{Tr}\left[\bm{K}_{YD}\bm{A}\bm{Z}\right]}{\partial\bm{K}_{YD}}\right)^{\top}\cdot\frac{\partial\bm{K}_{YD}}{\partial\bm{d_{ij}}}\right]\\\
&=&\operatorname{Tr}\left[\bm{A}\bm{Z}\cdot\frac{\partial\bm{K}_{YD}}{\partial\bm{d_{ij}}}\right]\end{array}$

and

$\begin{array}[]{rcl}\frac{\partial\operatorname{Tr}\left[\bm{Z}^{\top}\bm{A}^{\top}\bm{K}_{DD}\bm{A}\bm{Z}\right]}{\partial\bm{d_{ij}}}&=&\operatorname{Tr}\left[\left(\frac{\partial\operatorname{Tr}\left[\bm{Z}^{\top}\bm{A}^{\top}\bm{K}_{DD}\bm{A}\bm{Z}\right]}{\partial\bm{K}_{DD}}\right)^{\top}\cdot\frac{\partial\bm{K}_{DD}}{\partial\bm{d_{ij}}}\right]\\\
&=&\operatorname{Tr}\left[\bm{A}\bm{Z}\bm{Z}^{\top}\bm{A}^{\top}\cdot\frac{\partial\bm{K}_{DD}}{\partial\bm{d_{ij}}}\right].\end{array}$

The two partial derivatives of kernel matrices with respect to the current
atom are sparse matrices with non-zero columns or rows only where their index
is equal to the index of the current atom. The two matrices are computed as
follows:

$\frac{\partial\bm{K}_{YD}}{\partial\bm{d}_{ij}}=\left[\begin{array}[]{ccccc}0&\cdots&\frac{\partial
k\left(\bm{y}_{1},\bm{d}_{j}\right)}{\partial\bm{d}_{ij}}&\cdots&0\\\
\vdots&&\frac{\partial
k\left(\bm{y}_{2},\bm{d}_{j}\right)}{\partial\bm{d}_{ij}}&&\vdots\\\
\vdots&&\vdots&&\vdots\\\ 0&\cdots&\frac{\partial
k\left(\bm{y}_{N}\bm{d}_{j}\right)}{\partial\bm{d}_{ij}}&\cdots&0\end{array}\right]$
(6)

and

$\frac{\partial\bm{K}_{DD}}{\partial\bm{d}_{ij}}=\left[\begin{array}[]{ccccc}0&\cdots&\frac{\partial
k\left(\bm{d}_{1},\bm{d}_{j}\right)}{\partial\bm{d}_{ij}}&\cdots&0\\\
0&\cdots&\frac{\partial
k\left(\bm{d}_{2},\bm{d}_{j}\right)}{\partial\bm{d}_{ij}}&\cdots&0\\\
\vdots&&\vdots&&\vdots\\\ \frac{\partial
k\left(\bm{d}_{j},\bm{d}_{1}\right)}{\partial\bm{d}_{ij}}&\cdots&\frac{\partial
k\left(\bm{d}_{j},\bm{d}_{j}\right)}{\partial\bm{d}_{ij}}&\cdots&\frac{\partial
k\left(\bm{d}_{j},\bm{d}_{m}\right)}{\partial\bm{d}_{ij}}\\\
\vdots&&\vdots&&\vdots\\\ 0&\cdots&\frac{\partial
k\left(\bm{d}_{n},\bm{d}_{j}\right)}{\partial\bm{d}_{ij}}&\cdots&0\\\
\end{array}\right].$ (7)

Depending on the kernel function of interest, the partial derivatives are
computed via

$\frac{\partial
k\left(\bm{x},\bm{y}\right)}{\partial\bm{x}}=-\exp{\frac{-\|\bm{x}-\bm{y}\|^{2}_{2}}{2\sigma^{2}}}\frac{\left(\bm{x}-\bm{y}\right)}{\sigma^{2}}$

for the radial basis function kernel
$k(\bm{x},\bm{y})=\exp(\frac{-\left\|\bm{x}-\bm{y}\right\|_{2}^{2}}{\sigma^{2}})$,
and

$\frac{\partial
k\left(\bm{x},\bm{y}\right)}{\partial\bm{x}}=\beta\left(\bm{x}^{\top}\bm{y}+\alpha\right)^{\beta-1}\bm{y}$

for the polynomial kernel
$k(\bm{x},\bm{y})=(\bm{x}^{\top}\bm{y}+\alpha)^{\beta}$.

Since explicit solutions to the optimization problem on $\bm{d}_{j}$ do not
seem to exist, we update each column $\bm{d}_{j}$ through gradient descent
procedure applied on each element. At gradient descent iteration $\ell$, the
next element value is computed via

$\bm{d}_{ij}^{(\ell+1)}=\bm{d}_{ij}^{(\ell)}-\gamma\bm{g}_{ij}^{(\ell)},$

where $\gamma\in\mathbb{R}_{+}$ represents the chosen learning rate and
$\bm{g}_{ij}^{(\ell)}$ is the gradient

$\bm{g}_{ij}^{(\ell)}=\operatorname{Tr}\left[\bm{A}\bm{Z}\left(\bm{Z}^{\top}\bm{A}^{\top}\cdot\frac{\partial\bm{K}_{DD}}{\partial\bm{d_{ij}}}-2\cdot\frac{\partial\bm{K}_{YD}}{\partial\bm{d_{ij}}}\right)\right],$

where the matrices (6) and (7) are computed for the dictionary at iteration
$\ell$. The update step of the ORKDL-D algorithm is obtained by introducing a
few steps of gradient descent as described above for updating the dictionary
$\bm{D}$, after the update of $\bm{A}$ described by Algorithm 1.

###  III-C Mixed Optimized Reduced Kernel Dictionary Learning

The previous subsection presents an update procedure for the kernel vectors
$\bm{D}$. A possible drawback is that direct optimization of the objective of
(3) may lead to a dictionary $\bm{D}$ with weaker representation power for the
entire set of signals, $\bm{Y}$. We present here a further improvement. To
keep the representation significance of $\bm{D}$, which is a natural property
to require in the context of kernel methods, we introduce the standard DL
objective (1) into the ORKDL-D problem, obtaining

$\begin{array}[]{cl}\displaystyle\min_{\bm{A},\bm{Z},\bm{D},\bm{X}}&\|\varphi(\bm{Y})-\varphi(\bm{D})\bm{A}\bm{Z}\|_{F}^{2}+\lambda\|\bm{Y}-\bm{D}\bm{X}\|_{F}^{2}\\\
\text{ s.t. }&\left\|\bm{z}_{\ell}\right\|_{0}\leq s_{z},\ell=1:N\\\
&\left\|\varphi(\bm{D})\bm{a}_{j}\right\|=1,j=1:n\\\
&\left\|\bm{x}_{\ell}\right\|_{0}\leq s_{x},\ell=1:N\\\
&\left\|\bm{d}_{j}\right\|=1,j=1:n,\end{array}$ (8)

where $\lambda\in\mathbb{R}_{+}$ is a penalty constant. Since a mixed
objective is proposed, we name this method Mixed Optimized Reduced Kernel
Dictionary Learning (MORKDL-D).

By following two directions of optimization, we update the nonlinear
representation while conserving the representation power of the linear space.
The current optimization problem is solved similarly with the previous
problems. The update of a column of $\bm{D}$ is slightly different. The
gradient of the current atom uses also the partial derivative of the linear
term and is

$\tilde{\bm{g}}_{ij}=\bm{g}_{ij}-\bm{e}_{ij}$

where $\bm{e}_{ij}$ represents the $i$th element of the vector
$(\bm{Y}-\bm{D}\bm{X})\cdot\bm{X}_{j}^{\top}$; here, $\bm{X}_{j}^{\top}$ is
the $j$th row of $\bm{X}$. On the other hand, the linear sparse representation
matrix $\bm{X}$ is computed with the OMP algorithm. As before, the nonlinear
sparse representation matrix $\bm{Z}$ is computed using the Kernel OMP
algorithm, since $\bm{Z}$ and $\bm{X}$ can be optimized independently when
$\bm{A}$ and $\bm{D}$ are fixed.

##  IV Experiments

In this section we present the main results obtained with the proposed
methods. We used three different data sets: Digits, MNIST and CIFAR-10. For
each experiment we trained a kernel dictionary for representing the whole data
set or a subset extracted by its corresponding label. For example, for the
Digits data set we used all the signals during the training procedure, while
for the MNIST and CIFAR-10 data sets we used signals specific to a selected
label. For the MNIST data set we used label $5$, while for the CIFAR-10 data
set we used the first label.

All the algorithms were implemented in Matlab and Python and were run on a
Desktop PC with Ubuntu 20.04 as operating system. The PC has a processor with
16 cores, with a base frequency of 2.90 GHz (Max Turbo Frequency 4.80 GHz),
and 80 GB RAM memory. During the experiments, we computed the objective
function of (3) (we report the values
$\|\varphi(\bm{Y})-\varphi(\bm{D})\bm{A}\bm{Z}\|_{F}/\sqrt{mN}$ of the error
per signal element) at each iteration and measured the overall execution time.
Note that for MORKDL-D we report the values of the same error, not of the
objective of (8). Of course, for KDL we compute the objective of (2), which
corresponds to the standard case $\bm{D}=\bm{Y}$. The error and execution time
were computed by the mean on ten different rounds. As kernel function we used
the radial basis function
$k(\bm{x},\bm{y})=\exp(\frac{-\left|\bm{x}-\bm{y}\right|_{2}^{2}}{\sigma^{2}})$.
The kernel parameters were chosen according to a grid search, and for all data
sets we chose $\sigma=10$. All the algorithms trained a kernel dictionary
$\bm{A}$ of size $p=20$, with sparsity level $s_{z}=4$, using $10$ iterations.
For the reduced methods we initially trained a dictionary $\bm{D}$ of size
$n=50$, having a sparsity level of $s_{x}=5$, with $10$ iterations of the AK-
SVD algorithm. The dictionary $\bm{D}$ dimensions ensures a reduced kernel by
having a smaller number of atoms compared to the number of data signals. For
example, the Digits data set contains $N=5000$ signals of size $m=784$, while
the MNIST data set consists of approximately $N=6000$ signals of the same
size, of the same class. From the CIFAR-10 data set we used $N=5000$ signals,
specific to the interest class, of size $m=1024$.

In the ORKDL-D and MORKDL-D algorithms, the dictionary $\bm{D}$ is updated
with three gradient descent iterations. For each of the three data sets we
chose a learning rate that ensures a smooth decrease of the objective
function. We took $\gamma=5\cdot 10^{-4}$ for the Digits and MNIST data sets
and $\gamma=6\cdot 10^{-4}$ for the CIFAR-10 data set. For all our experiments
we used $\lambda=1$.

For a better visualization of the results, we show the evolution of the
nonlinear error (objective function of (3)) during the training procedure for
each algorithm using the three data sets Digits (Figure 1), MNIST (Figure 2)
and CIFAR-10 (Figure 3). As we can see, the improvements are visible from
iteration two of training. More results are given in two tables, containing
the value of the error at the last iteration (Table I) and the execution time
(Table II). All the reduced methods achieve a smaller error compared to the
KDL method. At first glance we notice that very large kernel spaces are not
necessary to obtain good results. The introduction of a small dimensional
space is enough to obtain satisfactory results. For situations where an
improvement of the nonlinear representation space is needed, several
iterations of gradient descent can be run as needed. According to the results
obtained on the three data sets it can be seen that the biggest advantage of
the proposed methods is the reduction of the execution time. For example, the
non-optimized reduced form obtains an execution time at least six times
shorter than the KDL method. In the problems where we try to train the kernel
space, an additional time will be introduced (for dictionary $\bm{D}$ update),
but the execution times are still shorter than those of KDL. For reduced
kernel methods the execution time is reduced by approximately 20-45% with
respect to the KDL problem. All the implementations are available at
https://github.com/denisilie94/rkdl.

Figure 1: Digits - representation error per iteration

Figure 2: MNIST - representation error per iteration

Figure 3: CIFAR-10 - representation error per iteration TABLE I: Last error value Algorithm Data set | Digits | MNIST | CIFAR-10
---|---|---|---
KDL | 1.041 $\cdot 10^{-2}$ | 8.081 $\cdot 10^{-3}$ | 3.373 $\cdot 10^{-3}$
RKDL-D | 1.026 $\cdot 10^{-2}$ | 7.970 $\cdot 10^{-3}$ | 3.285 $\cdot 10^{-3}$
ORKDL-D | 1.022 $\cdot 10^{-2}$ | 7.859 $\cdot 10^{-3}$ | 3.250 $\cdot 10^{-3}$
MORKDL-D | 1.021 $\cdot 10^{-2}$ | 7.846 $\cdot 10^{-3}$ | 3.238 $\cdot 10^{-3}$
TABLE II: Execution time in seconds Algorithm Data set | Digits | MNIST | CIFAR-10
---|---|---|---
KDL | 61.7 | 92.5 | 65.1
RKDL-D | 9.2 | 10.9 | 9.6
ORKDL-D | 39.2 | 48.4 | 51.6
MORKDL-D | 38.9 | 49 | 51.4

##  V Conclusions

In this paper we have presented a new approach to the Kernel Dictionary
Learning problem by introducing a reduced kernel and thus obtaining a new
algorithm, named RKDL. This method is suitable for problems with large data
sets, where standard KDL requires large memory sizes and long execution times.
Moreover, we have demonstrated that most of the time a large representation
space for the kernel matrix is not always needed to obtain satisfactory
results. The reduced form of the kernel is enough to get similar results or
even better. Since the kernel matrix is smaller, the required execution time
is also much shorter.

The RKDL algorithm was presented under three different forms: standard RKDL-D,
optimized RKDL-D and mixed optimized RKDL-D, the latter two including
optimization of the kernel vectors. All of them obtain competitive results
with the standard KDL problem.

## References

  * [1] H. Van Nguyen, V.M. Patel, N.M. Nasrabadi, and R. Chellappa.  Design of non-linear kernel dictionaries for object recognition.  IEEE Transactions on Image Processing, 22(12):5123–5135, 2013.
  * [2] J.J. Thiagarajan, K.N. Ramamurthy, and A. Spanias.  Multiple kernel sparse representations for supervised and unsupervised learning.  IEEE Transactions on Image Processing, 23(7):2905–2915, 2014.
  * [3] A. Golts and M. Elad.  Linearized kernel dictionary learning.  IEEE Journal of Selected Topics in Signal Processing, 10(4):726–739, 2016.
  * [4] S. Kumar, M. Mohri, and A. Talwalkar.  Sampling methods for the Nyström method.  The Journal of Machine Learning Research, 13(1):981–1006, 2012\.
  * [5] A. Rahimi and B. Recht.  Random features for large-scale kernel machines.  Advances in neural information processing systems, 20, 2007.
  * [6] Y.J. Lee and O.L. Mangasarian.  RSVM: Reduced support vector machines.  In Proceedings of the 2001 SIAM International Conference on Data Mining, pages 1–17. SIAM, 2001.
  * [7] K.M. Lin and C.J. Lin.  A study on reduced support vector machines.  IEEE Transactions on Neural Networks, 14(6):1449–1459, 2003.
  * [8] G.M. Fung, O.L. Mangasarian, and A.J. Smola.  Minimal kernel classifiers.  Journal of Machine Learning Research, 3(Nov):303–321, 2002.
  * [9] W. Deng, Q. Zheng, and K. Zhang.  Reduced kernel extreme learning machine.  In Proceedings of the 8th international conference on computer recognition systems CORES 2013, pages 63–69. Springer, 2013.
  * [10] J. Hu and Y.P. Tan.  Nonlinear dictionary learning with application to image classification.  Pattern Recognition, 75:282–291, 2018.
  * [11] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [12] A. Krizhevsky, G. Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * [13] Y.C. Pati, R. Rezaiifar, and P.S. Krishnaprasad.  Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition.  In Proceedings of 27th Asilomar conference on signals, systems and computers, pages 40–44. IEEE, 1993.
  * [14] B. Dumitrescu and P. Irofti.  Dictionary learning algorithms and applications.  Springer, 2018.
  * [15] M. Aharon, M. Elad, and A. Bruckstein.  K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation.  IEEE Transactions on Signal Processing, 54(11):4311–4322, 2006\.
  * [16] R. Rubinstein, M. Zibulevsky, and M. Elad.  Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit.  Technical report, Computer Science Department, Technion, 2008.
  * [17] D.C. Ilie-Ablachim and B. Dumitrescu.  Anomaly detection with selective dictionary learning.  In International Conference on Control, Decision and Information Technologies. IEEE, 2022.