# Dynamic Memory Based Adaptive Optimization

Balázs Szegedy†, Domonkos Czifra†, Péter Kőrösi-Szabó†

†Alfréd Rényi Insititute of Mathematics, Budapest, Hungary
{szegedy, doma945<EMAIL_ADDRESS>
###### Abstract

Define an optimizer as having memory $k$ if it stores $k$ dynamically changing
vectors in the parameter space. Classical SGD has memory $0$, momentum SGD
optimizer has $1$ and Adam optimizer has $2$. We address the following
questions: How can optimizers make use of more memory units? What information
should be stored in them? How to use them for the learning steps? As an
approach to the last question, we introduce a general method called
”Retrospective Learning Law Correction” or shortly RLLC. This method is
designed to calculate a dynamically varying linear combination (called
learning law) of memory units, which themselves may evolve arbitrarily. We
demonstrate RLLC on optimizers whose memory units have linear update rules and
small memory ($\leq 4$ memory units). Our experiments show that in a variety
of standard problems, these optimizers outperform the above mentioned three
classical optimizers. We conclude that RLLC is a promising framework for
boosting the performance of known optimizers by adding more memory units and
by making them more adaptive.

##  1 Introduction

In this paper, we investigate optimizers that store $k$ vectors in the
parameter space $\mathbb{R}^{n}$ of a neural network or more generally in the
parameter space related to any optimization problem. We call such vectors
memory units and we measure the memory usage of an optimizer by the number
them.

The simplest example for an optimizer with memory is the momentum SGD
optimizer which stores a single vector $m$ (momentum vector) in the parameter
space $\mathbb{R}^{n}$. In each step, $m$ is updated according to the rule
$m\longleftarrow\beta m+\bigtriangledown_{\theta}f(\theta)$ where $f$ is the
objective function, $\theta\in\mathbb{R}^{n}$ is the parameter vector and
$0<\beta<1$ is a fixed real number. The vector $\theta$ is updated according
to the rule $\theta\longleftarrow\theta-cm$ where $c>0$ is the learning rate.

The Adam optimizer operates with two memory units. One of them is the momentum
vector and the other one is the momentum of the squares of the gradient
vectors. In contrast with the momentum optimizer, the Adam optimizer is not
linear in the gradient vectors. Neither the update rule of the memory units,
nor the way the memory units are used for the parameter update is linear.

The present paper has two independent contributions. The first contribution is
a novel and simple method that we call RLLC=Retrospective Learning Law
Correction. It is an update rule for a vector $L$ (called learning law) that
describes a natural way of using a set of dynamically changing memory units
for the update of the parameter vector $\theta$. More precisely,
$L\in\mathbb{R}^{k}$ contains the coefficients of a linear combination of the
$k$ memory units which is multiplied by a fixed learning rate $c_{1}$ and
substracted from $\theta$ as usual. In each step, before updating $\theta$ and
the memory units, we update $L$ by the formula $L\longleftarrow L+c_{2}M^{+}g$
where $M$ is the $n\times k$ matrix formed by the memory units, $M^{+}$ is the
Moore-Penrose inverse of $M$, $g$ is the newly received gradient and $c_{2}$
is the meta learning rate. The main idea behind the update rule for $L$ is
that the new gradient $g$ contains retrospective information on how the
algorithm could have performed better in the previous step. Thus it can be
used to compute a corrected version of $L$ which ”thinks more ahead”. Note
that our update rule of the learning law can also be regarded as a general
framework for associating a $k$-dimensional adaptive learning rate with an
arbitrary set of $k$ evolving memory units.

As the second main contribution, we examine optimizers in which memory units
are updated by fixed linear rules. More precisely, in each step each memory
unit is updated to a linear combination of the memory units and the new
arriving gradient. The parameter vector is updated by a (possibly changing)
liner combination (given by the learning law $L$) of the memory units. Such
optimizers are interesting even if the learning law is fixed. They include
SGD, momentum SGD and Nesterov Accelerated Gradient (NAG) Nesterov (2012).
Thus, the linear framework provides a useful generalization of these famous
optimizers and enables a dynamically changing continuous interpolation between
them. The RLLC method turns out to be ideal for this. Our experiments show
that linear memory combined with RLLC leads to powerful optimizers. The case
of memory $1$ is already interesting. A memory $1$ linear optimizer stores a
momentum vector. Applying RLLC in this trivial setting yields a variant of the
momentum SGD optimizer enhanced with a new type of adaptive learning rate. As
the number of memory units increases, the mathematics becomes more complex,
presenting a field of study that is interesting in its own right. We present
some of the fundamental properties of linearly updated memory units. In
particular, we prove a version of basis independence for RLLC combined with
linear memory which allows us to apply basis transformations to the update
rules without changing the optimization process. This together with a variant
of the Jordan normal form over the field $\mathbb{R}$ helps to convert these
optimizers into a canonical form in which each memory unit is associated with
a so-called Jordan block. A Jordan block of size $1$ corresponds to a single
memory unit (denoted by $M(\beta)$) storing a momentum vector of the gradients
with parameter $\beta$. A Jordan block of size $2$ either corresponds to a
pair of memory units (denoted by $CM(\beta),\beta\in\mathbb{C}$) namely the
real and imaginary parts of a momentum vector with complex parameter or to a
pair of memory units $m_{1},m_{2}$ (denoted by $M_{2}(\beta)$) where $m_{1}$
is a momentum vector of the gradient and $m_{2}$ is a momentum vector of
$m_{1}$, both with parameter $\beta$. In general, there are two infinite
families of Jordan blocks giving rise to $k$-tuples or $2k$-tuples of memory
units denoted by $M_{k}(\beta)$ and $CM_{k}(\gamma)$. These are the
fundamental building blocks of linearly updated memory. We denote the natural
operation by $\oplus$ which combines these building blocks into larger memory
by the union of the corresponding memory units. By slightly abusing the
notation we often identify memory update rules with optimizers where learning
is given by the RLLC method. For example, $M(\beta)$ also denotes the memory
$1$ optimizer with memory unit $M(\beta)$ (a momentum vector) and with RLLC.
Notice that the $M(\beta)$ optimizer is a close relative of momentum SGD but
it is not equivalent with it.

In our experiments, we identified a number of interesting simple settings
involving few (at most $4)$ memory units. These include the types of
optimizers $M(\beta)$, $M(\beta)\oplus M(0)$, $M(\beta_{1})\oplus
M(\beta_{2})\oplus M(\beta_{3})$, $M_{2}(\beta)$, $M_{3}(\beta)$ and
$M(\beta){\oplus}M(-\beta){\oplus}CM(\beta i)$. We observed that these
optimizers often surpassed the performance of three commonly used optimizers
across a variety of tasks even without carefully optimizing the parameters
$\beta_{i}$. Notice that $M(\beta)\oplus M(0)$ is an adaptively changing
linear combination of SGD, momentum SGD and NAG. Thus, it adaptively
interpolates between three well known optimizers. Remarkably, it demonstrated
competitive or even superior performance compared to the Adam optimizer in
many tasks, which also uses two memory units.

This paper primarily aims not to challenge all existing optimizers in the
field, but rather to introduce a novel mathematical concept that could spark
further research. The experimental results presented here should be
interpreted as an illustration of the potential of our approach. We posit that
the implications of the RLLC method extend beyond mere enhancements to current
optimization techniques, suggesting broader applications and insights in the
realm of optimization and machine learning.

##  2 Related Work

In the field of deep neural networks, the most commonly used optimizers are
still hand-designed, such as momentum SGD, RMSProp, or Adam Kingma and Ba
(2014). Tremendous attempts were made to improve their performance, or exceed
them with different approaches. Finding the best optimizer in general is not
possible, see Wolpert and Macready (1997). However by restricting the tasks,
and hence, finding a better optimizer can be regarded as a learning process.
This learning can focus either within a single task, or across multiple tasks.

For the former, typical approaches include adaptive learning based methods
(Wang et al., 2022; Keskar and Socher, 2017b) and extensions of hand-crafted
optimizers (e.g.: by replacing some components with RNN) (Li and Malik, 2017a,
b; Chen et al., 2022). However these methods aim to speed up single task
convergence, but adaptive learning may have some drawbacks, especially on NLP
tasks, and certain synthetic problems Wilson et al. (2018).

The later one, namely learned optimization, trains a function (typically an
auxiliary neural network) to optimize the model that solves the original
tasks. Some of them employ traditional gradient-descent based optimizers
(Andrychowicz et al., 2016; Metz et al., 2022b; Bengio et al., 1991;
Wichrowska et al., 2017), other methods utilize evolutionary algorithms
(Bengio et al., 1991; Metz et al., 2020), or reinforcement learning (Li and
Malik, 2017c; Bello et al., 2017). These methods aim to discover an optimizer
via some form of learning that enhance performance over a diverse set of
tasks, surpassing the effectiveness of commonly used, hand-designed
optimizers.

Our approach aligns more with the adaptive learning based methods. However,
rather than focusing on performance enhancement for specific group of tasks,
our objective is to establish a comprehensive mathematical framework, which
supports the combination of many existing optimizers, and enables the
exploration of new optimization algorithms. While our goal is different, our
work shows some similarities with learned-optimization as well. Training an
optimizer involves backpropagating through many unrolled inner optimization
steps, a process that is not only computationally intensive but also
susceptible to instability, as highlighted by Metz et al. (2020). Within this
context, our RLLC method could be interpreted as a kind of differential-free
optimization technique. In this paper, we focus solely looking one step ahead,
but one could easily extend it beyond that.

##  3 Retrospective Learning Law Correction

Functional approach to optimizers: The RLLC method is presented through an
abstract mathematical framework for optimizers. This framework is somewhat
specialized, yet it maintains sufficient generality to encompass a range of
interesting optimizers. We think of optimizers as entities with an evolving
internal state that updates at each step based on newly received gradients.
Additionally, the optimizer calculates a parameter update vector relevant to
the optimization process. A functional description of such an optimizer is
given in the following definition.

######  Definition 3.1.

An optimizer for $n$ parameters is a pair of functions of the form
$F:\mathcal{S}\times\mathbb{R}^{n}\to\mathcal{S}$ and
$G:\mathcal{S}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ where $\mathcal{S}$ is
the set of possible internal states, $F$ is the state update function and $G$
is the parameter update function.

To translate optimizers into an actual optimization process we choose an
initial internal state $S_{0}\in\mathcal{S}$ and an initial parameter vector
$\theta_{0}\in\mathbb{R}^{n}$. Then we iterate

$S_{t}:=F(S_{t-1},g_{t}),~{}~{}\theta_{t}:=\theta_{t-1}+G(S_{t},g_{t})$

where $g_{t}$ is a gradient vector received by the optimizer in the $t$-th
step. To illustrate this formalism, assume that the optimizer is given by
$\mathcal{S}=\mathbb{R}^{n}$, $F(v,w)=\beta v+w,G(v,w)=-cv$. In this case, we
obtain the momentum SGD with learning rate $c$ and decay parameter $\beta$.

Optimizers with memory and RLLC: We will think of memory $k$ optimizers in a
way that the internal state space is of the form $\mathbb{R}^{n\times
k}\times\mathcal{H}$ where the columns of matrices in $\mathbb{R}^{n\times k}$
represent $k$ vectors in the parameter space $\mathbb{R}^{n}$ and
$\mathcal{H}$ will be called the space of hidden states. We typically assume
that $n$ is a large number and that the hidden states are described by much
fewer than $n$ parameters. A memory update rule for $k$ memory units is a
function of the form $U:(\mathbb{R}^{n\times
k}\times\mathcal{H})\times\mathbb{R}^{n}\to\mathbb{R}^{n\times
k}\times\mathcal{H}$ where the external $\mathbb{R}^{n}$ component represents
new arriving gradients. Such a function does not yet determine an optimizer.
The RLLC method is designed to turn memory update rules into optimizers by
extending their state space with a vector called learning law and introducing
a natural parameter update function. We give two different descriptions of
RLLC. The first one is a functional description which is more convenient for
proofs.

We will need the so-called Moore-Penrose inverse which is defined for an
arbitrary matrix $A\in\mathbb{R}^{n\times k}$ and is denoted by $A^{+}$. Note
that if $A$ has rank $k$ (which means that $A$ is non-degenerate if $n\geq k$)
then $A^{+}=(A^{T}A)^{-1}A^{T}$.

######  Definition 3.2 (RLLC functional form).

Let $U$ be a memory update rule as above. Then the corresponding RLLC
optimizer with learning rates $c_{1},c_{2}$ is given as follows. The state
space is $\mathcal{S}:=\mathbb{R}^{n\times
k}\times\mathcal{H}\times\mathbb{R}^{k}$ where the extra component
$\mathbb{R}^{k}$ is called the learning law. The functions $F,G$ are given in
the following way. Assume that $M\in\mathbb{R}^{n\times
k},H\in\mathcal{H},L\in\mathbb{R}^{k},g\in\mathbb{R}^{n}$. Then

$F(M,H,L,g):=(U_{1}(M,H),U_{2}(M,H),L+c_{2}M^{+}g)$ $G(M,H,L,g):=-c_{1}ML.$

######  Remark 3.3.

Vectors in $\mathbb{R}^{m}$ are considered to be column vectors. This means
that they are treated as $m\times 1$ matrices in calculations.

The second, less abstract approach to RLLC describes the optimization process
directly in a more conventional way.

Require:

  * •

$\theta_{0}$ initial parameter vector

  * •

$f(\theta)$: stochastic objective function with parameters
$\theta\in\mathbb{R}^{n}$

  * •

Two learning rates $c_{1},c_{2}>0$

  * •

Stability parameter $\epsilon$ for relaxed Penrose inverse

  * •

$M_{0}\in\mathbb{R}^{n\times k}$: initial memory units

  * •

$L_{0}\in\mathbb{R}^{k}$: initial learning law

  * •

$H_{0}$: initial hidden state

$t:=0$: initialize time step

while: $\theta_{t}$ not converged do:

$t\longleftarrow t+1$

$g_{t}:=\bigtriangledown_{\theta}f_{t}(\theta_{t-1})$ (Get gradients w.r.t.
stochastic objective at timestep t) Get New Gradient

$L_{t}:=L_{t-1}+c_{2}M_{t-1}^{+}g_{t}$ RLLC Step

$(M_{t},H_{t}):=U(M_{t-1},H_{t-1},g_{t})$ Memory Update

$\theta_{t}:=\theta_{t-1}-c_{1}M_{t}L_{t}$ Learning Step

Explanation of RLLC and remarks:

The idea behind the learning rule update is that with the arrival of the new
gradient $g_{t}$ the optimizer gains new (retrospective) information on how it
could have done better in the previous learning step. Notice that the vector
$M^{+}g_{t}$ is the coefficient vector of the orthogonal projection of $g_{t}$
to the space spanned by the memory units when written as a linear combination
of the memory units. This means that, if performed with the new law, the
outcome of the parameter update in step $t-1$ would have been
$\theta_{t}-c_{1}c_{2}p_{t}$ instead of $\theta_{t}$ where $p_{t}$ is the
orthogonal projection of $g_{t}$ to the space spanned by the memory units in
the $t-1$-th step. Notice that $(p_{t},g_{t})=(p_{t},p_{t})\geq 0$ and thus
the change $-c_{1}c_{2}p_{t}$ points in a direction which improves the
objective function.

The above heuristics does not take it into account that the objective function
$f_{t}$ is also changing. This fact indicates that our update rule is more
justified if the second learning rate $c_{2}$ is small and thus random effects
have time to average out leaving only useful directions in the update. Also
notice that the algorithm does not ”go back in time” to perform the improved
learning step. Instead it applies the updated learning law with the updated
memory units. This shows that the efficiency of the RLLC method depends on a
type of consistency property. Roughly speaking it assumes that the notion of a
”good learning law” does not change too much in time and so improvements of
the past give improvements of the future. For this reason the choice of the
memory update rule is a crucial issue which is one of the main topics of the
second part of this paper.

######  Remark 3.4.

The performance of the RLLC optimizer is dependent on the initialization of
the learning law at the beginning. In practice it is not initialized to be
$0$.

######  Remark 3.5.

To avoid numerical instability, in practice we use a relaxed version of
Penrose inverse which has a parameter $\epsilon$ set to a small number.

Linear invariance of RLLC: We close this chapter with a useful linear
invariance property of the RLLC method. We will need the next two definitions.

######  Definition 3.6.

Let $U$ be a memory update rule as above and let $Q\in\mathbb{R}^{k\times k}$
be an invertible matrix. We define the new memory update rule $U^{Q}$ in the
following way: Let $M\in\mathbb{R}^{n\times k}$, $H\in\mathcal{H}$ and
$U(MQ^{-1},H,g)=(M_{2},H_{2})$. Then

$U^{Q}(M,H,g):=(M_{2}Q,H_{2}).$

######  Definition 3.7.

Two optimizers given by $(F,G)$ and $(F^{\prime},G^{\prime})$ with state
spaces $\mathcal{S}$ and $\mathcal{S}^{\prime}$ are called equivalent if there
is a bijection $\phi:\mathcal{S}\to\mathcal{S}^{\prime}$ (called an
isomorphism) such that $\phi(F(S,v))=F^{\prime}(\phi(S),v)$ and
$G(S,v)=G^{\prime}(\phi(S),v)$. A partial isomorphism is a bijection between a
subset of $\mathcal{S}$ and a subset of $\mathcal{S}^{\prime}$ having the same
property. If there is such a function we say that the two optimizers are
partially equivalent on these two subsets. In particular, if two optimizers
with memory $k$ states are partially equivalent on states with rank $k$ memory
matrices then we call them essentially equivalent.

It is easy to see that if two optimizers are equivalent then they define the
same optimization process if their initialization of internal states is
isomorphic. If two optimizers are partially equivalent with partial
isomorphism $\phi$ then the optimization processes are identical as long as
they operate on states in the domain (and image) of $\phi$.

######  Lemma 3.8 (Linear invariance of RLLC).

Let $U$ be a memory update rule as above and $Q\in\mathbb{R}^{k\times k}$ be
an arbitrary matrix. Then the RLLC optimizer corresponding to $U$ is
essentially equivalent to the RLLC optimizer corresponding to $U^{Q}$.

###### Proof.

We claim that the function $\phi(M,H,L):=(MQ,H,Q^{-1}L)$ is a partial
isomorphism on states with rank $k$ memory matrices. This follows trivially
from formulas in definition 3.2 and the fact that $(MQ)^{+}=Q^{-1}M^{+}$ holds
if ${\rm rank}(M)=k$. ∎

##  4 Linear memory updates

Throughout this chapter we investigate linear memory update rules with no
hidden states. Such an update rule is given by

$U(M,g):=MB+ga^{T}$ (1)

where $M\in\mathbb{R}^{n\times k}$ is the memory matrix, $g\in\mathbb{R}^{n}$
is a new gradient and $a\in\mathbb{R}^{k},B\in\mathbb{R}^{k\times k}$ are
fixed parameters of the update rule. In an optimization process this means
that the memory unit $m_{i}$ represented by the $i$-th column of $M$ is
updated to $a_{i}g+\sum_{j=1}^{k}B_{j,i}m_{i}$ when the new gradient $g$ is
received.

Linear memory optimizers with fixed learning law: Linearly updated memory
units are interesting independently of the RLLC method. We can directly obtain
powerful optimizers by using a fixed hand designed learning law
$L\in\mathbb{R}^{k}$. This type of optimizer, denoted by $\mathcal{L}(B,a,L)$
works by the equations:

$M_{t}=M_{t-1}B+g_{t}a^{T}~{},~{}\theta_{t}=\theta_{t-1}-M_{t}L.$

If $k=1$ then $B,a,L$ are single real numbers. The corresponding optimizer is
a momentum SGD optimizer with decay parameter $B$ and learning rate $aL$.
Another, important setting is described in the following lemma (for proof see
Appendix A).

######  Lemma 4.1.

Let

$B=\begin{pmatrix}\beta&0\\\ 0&0\end{pmatrix},~{}a=\begin{pmatrix}1\\\
1\end{pmatrix},~{}L=\begin{pmatrix}c\beta\\\ c\end{pmatrix}$

Then the corresponding optimizer $\mathcal{L}(B,a,L)$ is Nestorv Accelerated
Gradient with decay parameter $\beta$ and learning rate $c$.

Abstract rule of a memory unit: There is a useful observation which sheds more
light on what information linear memory stores if this memory update is
iterated in an optimization process started with initial value $0^{n\times k}$
for $M$. By induction we have

$U(...U(U(0^{n\times
k},g_{1}),g_{2}),...,g_{t})=\sum_{i=1}^{t}g_{t}a^{T}B^{t-i}.$

We obtain that at time $t$ (after $t$ iteration of the update rule) the value
of the $i$-th memory unit is given by

$m_{j}=\sum_{i=0}^{t-1}g_{t-i}(a^{T}B^{i})_{j}$ (2)

where $(a^{T}B^{i})_{j}$ denotes the $j$-th coordinate of the row vector
$a^{T}B^{i}$. If we regard gradients with index $0$ or negative index as $0$
then the sum can be taken from $0$ to infinity. Informally speaking, this
means that $m_{i}$ is a fixed (time independent) linear combination of
previous gradients going backwards in time. This linear combination is
represented by the infinite sequence $\\{(a^{T}B^{i})_{j}\\}_{i=0}^{\infty}$
for the $j$-th memory unit. We say that this infinite sequence is the abstract
rule of the memory unit. To guarantee that older gradients are taken with
decaying weight in (2) we need to assume that the spectral norm of $B$ is
smaller than $1$.

Real momentum: Memory updates in the case $k=1$ are determined by two numbers:
$\beta=B_{1,1}$ and $\alpha=a_{1}$. The update rule of the single memory unit
$m$ in the $t$-th step is $m\longleftarrow\alpha g_{t}+\beta m$. (This is
essentially the update rule of a momentum vector.) It follow from our formula
that the abstract rule in this case is given by the geometric sequence
$\beta^{i}\alpha$. In particular, in the $t$-th step we have that
$m=\sum_{i=0}^{t}\alpha\beta^{i}g_{t-i}$.

Complex momentum: The case

$B=\begin{pmatrix}\alpha&-\beta\\\ \beta&\alpha\end{pmatrix}$ (3)

has a distinguished role because such matrices represent complex numbers
$\gamma=\alpha+\beta i$. In this special case the two memory units can be
interpreted as the real and the complex parts of a single complex valued
memory unit which describes a momentum vector with complex parameter $\gamma$.
More precisely $m_{1}$ and $m_{2}$ are the real and the complex parts of a
memory unit $m\in\mathbb{C}^{n}$ which is updated according to
$m\longleftarrow g_{t}+\gamma m$.

Jordan block of size $2$: Another interesting example for $k=2$ is given by

$B=\begin{pmatrix}\alpha&1\\\ 0&\alpha\end{pmatrix}$

which is the so called Jordan block of size $2$ with eigenvalue $\alpha$. The
first one of the two memory units is the momentum vector of the gradients with
parameter $\alpha$. However the second memory unit stores something new. It is
the momentum vector of the first memory unit with parameter $a$. One can show
that the abstract rule corresponding to this memory unit is given by the
infinite sequence $0,1,2\alpha,3\alpha^{2},4\alpha^{3},\dots$.

Propagators and their unions: In general, we call $k$ memory units
$m_{1},m_{2},\dots,m_{k}$ connected by a joint linear update rule a propagator
of dimension $k$. Recall that such an object is described by a matrix
$B\in\mathbb{R}^{k\times k}$ and a vector $a\in\mathbb{R}^{k}$. In the
previous two examples each real number $\beta$ is associated with a propagator
denoted by $M(\beta)$ of dimension $1$ while complex numbers $\gamma$ are
associated with a propagator denoted by $CM(\gamma)$ of dimension $2$. We call
such propagators momentum propagators. It will be important for us that there
is a simple operation on propagators that we call union and denote by
$\oplus$. This is simply just taking the union of the corresponding memory
units together and updating them independently. From a linear algebraic point
of view, the matrix $B$ corresponding to the union of propagators is a block
diagonal matrix whose blocks contain the matrices of the individual
propagators. The vector $a$ corresponding to the union is the concatenation of
the vectors of the propagators. Unions of momentum propagators will be called
multi momentum propagators.

##  5 Optimizers with linear memory and RLLC

In this chapter we discuss the basic properties of optimizers which combine
linear memory and the RLLC method. We use the term LM-RLLC optimizers for
them. Based on definition 3.2 and formula (1) one can produce the LM-RLLC
optimizer $\mathcal{F}(B,a,c_{1},c_{2})$ with hyperparameters
$B,a,c_{1},c_{2}$. The corresponding update functions are given by

$F(M,L,g)=(MB+ga^{T},L+c_{2}M^{+}g)$ $G(M,L,g)=-c_{1}ML$

For the sake of completeness we describe the recursive optimization process.

######  Definition 5.1.

(LM-RLLC optimization process) Let us fix the hyperparameters
$B\in\mathbb{R}^{k\times
k},a\in\mathbb{R}^{k},(c_{1},c_{2})\in\mathbb{R}^{2}.$ Then the LM-RLLC
optimizer with these hyperparameters is given by the equations

$L_{t}=L_{t-1}+c_{2}M_{t-1}^{+}g_{t}$ $M_{t}=M_{t-1}B+g_{t}a^{T}$
$\theta_{t}=\theta_{t-1}-c_{1}M_{t}L$

where $M_{0}$ is the $0$ matrix in $\mathbb{R}^{n\times k}$ and
$L_{0}\in\mathbb{R}^{k}$ is a suitable (typically non $0$) vector.

Deeper mathematical analysis reveals that LM-RLLC optimizers can be
transformed into a simpler, canonical form if we look at them up to
equivalence. The key observation is a ”basis independence” property of LM-RLLC
optimizer functions.

######  Theorem 5.2 (Basis independence of LM-RLLC optimizers).

Let $k\in\mathbb{N},a\in\mathbb{R}^{k},B\in\mathbb{R}^{k\times
k},(c_{1},c_{2})\in\mathbb{R}^{2}$ and let $Q\in\mathbb{R}^{k\times k}$ be an
invertible matrix. Then $\mathcal{F}(B,a,c_{1},c_{2})$ is essentially
equivalent to

$\mathcal{F}(Q^{-1}BQ,Qa,c_{1},c_{2}).$

###### Proof.

The optimizer $\mathcal{F}(B,a,c_{1},c_{2})$ is obtained from the linear
memory update rule $U(M,g)=MB+ga^{T}$ with RLLC. Notice that $U^{Q}$ (in the
sense of definition 3.6) is given by $U^{G}(M,g)=M(Q^{-1}BQ)+g(Qa)^{T}$. Then
lemma 3.8 finishes the proof. ∎

Real Jordan normal form: Theorem 5.2 together with a variant of the Jordan
decomposition theorem implies that we can transform LM-RLLC optimizers into a
very special form without changing the optimization process. The original form
of the Jordan decomposition theorem says that if $B\in\mathbb{C}^{k\times k}$
is an arbitrary complex matrix then there is an invertible matrix
$Q\in\mathbb{C}^{k\times k}$ such that $Q^{-1}BQ$ has a block diagonal form
with each block being a so-called Jordan block. A Jordan block
$J_{m}(\lambda)$ is a matrix of size $m\times m$ with $\lambda\in\mathbb{C}$
in the diagonal, $1$ above the diagonal and $0$ everywhere else. For example

$J_{3}(\lambda)=\begin{pmatrix}\lambda&1&0\\\ 0&\lambda&1\\\ 0&0&\lambda\\\
\end{pmatrix}$

There is a similar, although somewhat more complicated statement (called real
Jordan normal form) if $B$ and $Q$ are required to be real matrices. In this
case there are two types of blocks $J_{m}(\lambda)$ with
$\lambda\in\mathbb{R}$ and $CJ_{m}(\alpha+\beta i)$ with
$\alpha,\beta\in\mathbb{R}$. The second type of block has size $2m\times 2m$
and it ”imitates” complex Jordan blocks with real matrices. This matrix is
very similar to $J_{m}(\lambda)$ with the main difference being that each
entry is replaced by a $2\times 2$ matrix. The $0$’s and $1^{\prime}s$ are
replaced by $0$ matrices and identity matrices. The $\lambda$ entries are
replaced by the matrix in equation $(\ref{comprep})$ which represents
$\alpha+\beta i$ by a real matrix. For example

$CJ_{2}(\alpha+\beta i)=\begin{pmatrix}\alpha&-\beta&1&0\\\
\beta&\alpha&0&1\\\ 0&0&\alpha&-\beta\\\ 0&0&\beta&\alpha\end{pmatrix}$

Propagators of Jordan type: Racall that an LM-RLLC optimizer is given by
$B\in\mathbb{R}^{k\times k},a\in\mathbb{R}^{k}$ and two learning rates. By
transforming the matrix $B$ to its real Jordan normal form with a basis
transformation given by $Q\in\mathbb{R}^{k\times k}$, we can divide the memory
units into groups belonging to single blocks of type $J_{m}(\lambda)$ or
$CJ_{m}(\alpha+\beta i)$. The block diagonal form of $Q^{-1}BQ$ guarantees
that these groups do not interact with each other in memory updates and thus
we can treat them as separate propagators. Recall that in this basis
transformation considered in theorem 5.2 the vector $a$ transforms into $Qa$.
By applying a statement which is slightly stronger then the Jordan
decomposition theorem we can also guarantee that the part of $a$ in each block
contains at most one coordinate with $1$ and the rest is $0$. We can also
assume that this coordinate is the first one otherwise there are trivial
memory units which store $0$ in each step. By summarizing all of this we
obtain propagators of very special type. Let $e_{m}\in\mathbb{R}^{m}$ denote
the vector with $1$ in the first coordinate and $0$ in the rest. We denote the
propagator corresponding to the pair $(J_{m}(\lambda),e_{m})$ by
$M_{m}(\lambda)$ and the propagator corresponding to $(CJ_{m}(\alpha+\beta
i),e_{2m})$ by $CM_{m}(\alpha+\beta i)$. We call such propagators Jordan block
propagators. If $m=1$ then we omit the index and simply write $M(\lambda)$ and
$CM(\alpha+\beta i)$. We obtain the next theorem.

######  Theorem 5.3 (Normal forms of LM-RLLC optimizers).

Every LM-RLLC optimizer is essentially equivalent with another LM-RLLC
optimizer where the memory update is of the form $P_{1}\oplus
P_{2}\oplus\dots\oplus P_{r}$ where each $P_{i}$ is a Jordan block propagator.

By slightly abusing the notation we will also use the formula $P_{1}\oplus
P_{2}\oplus\dots\oplus P_{r}$ for the optimizer itself. For example
$M(0.9)\oplus M_{2}(0.6)\oplus CM_{2}(0.3+0.2i)$ stands for a memory $7$
optimizer where the memory units are grouped and updated according to the
propagators $M(0.9)$,$M_{2}(0.6)$ and $CM_{2}(0.3+0.2i)$.

##  6 Experiments

For our experiments, we used the Learned Optimization framework Metz et al.
(2022a) as a starting point. The framework offers pre-trained and hyper
parameter optimized optimizers. We compare our results with the most widely
used optimizers as baseline: Adam, SGD, and SGD with momentum. We compare test
loss and classification accuracy on MNISTDeng (2012), Fashin-MNISTXiao et al.
(2017b), and CIFAR-10Krizhevsky (2009) datasets. We experimented with dense,
convolutional and residual neural networks. The source code of our work is
available
publicly111https://anonymous.4open.science/r/easytrace-843E/README.md. See
implementation details in Appendix C.

| MNIST |  | Fashion-MNIST |  | CIFAR-10 |
---|---|---|---|---|---|---
| MLP | Conv |  | MLP | Conv |  | MLP | Conv | ResNet-20 |
$SGD$ | 0.0882 | 0.0331 |  | 0.3426 | 0.3673 |  | 1.4084 | 0.8359 | 0.6010 | Loss
98.16 | 98.56 |  | 88.65 | 86.97 |  | 52.18 | 71.37 | 80.93 | Acc
$MomentumSGD$ | 0.0856 | 0.0324 |  | 0.3476 | 0.2732 |  | 1.4108 | 0.7850 | 0.5757 | Loss
98.22 | 98.97 |  | 88.67 | 90.78 |  | 51.99 | 73.17 | 81.36 | Acc
$Adam$ | 0.0758 | 0.0304 |  | 0.3407 | 0.2704 |  | 1.3858 | 0.7920 | 0.5857 | Loss
97.83 | 98.99 |  | 88.78 | 90.76 |  | 52.43 | 73.70 | 81.48 | Acc
$M(0.9)$ | 0.0844 | 0.0310 |  | 0.3408 | 0.2661 |  | 1.4021 | 0.8030 | 0.5301 | Loss
98.21 | 99.03 |  | 88.64 | 90.96 |  | 52.13 | 73.95 | 83.08 | Acc
$M(0.9){\oplus}M(0.0)$ | 0.0888 | 0.0323 |  | 0.3475 | 0.2678 |  | 1.3973 | 0.7977 | 0.5353 | Loss
98.26 | 98.95 |  | 88.82 | 90.98 |  | 51.71 | 74.11 | 83.38 | Acc
$M(0.9){\oplus}M(0.8){\oplus}M(0.7)$ | 0.0829 | 0.0343 |  | 0.3359 | 0.2563 |  | 1.4142 | 0.7734 | 0.5268 | Loss
98.23 | 98.95 |  | 88.67 | 91.11 |  | 51.55 | 75.42 | 83.19 | Acc
$M_{2}(0.6)$ | 0.0801 | 0.0800 |  | 0.3220 | 0.4488 |  | 1.3444 | 1.0404 | 0.5811 | Loss
98.17 | 97.60 |  | 88.75 | 84.31 |  | 53.42 | 63.59 | 81.35 | Acc
$M(0.9){\oplus}M_{2}(0.6)$ | 0.0861 | 0.0319 |  | 0.3536 | 0.2636 |  | 1.4155 | 0.7602 | 0.5354 | Loss
98.21 | 98.99 |  | 89.03 | 90.95 |  | 52.08 | 75.87 | 82.84 | Acc
$M(0.9){\oplus}M(0.0){\oplus}M_{2}(0.6)$ | 0.0877 | 0.0287 |  | 0.3498 | 0.2596 |  | 1.4028 | 0.7216 | 0.5393 | Loss
98.23 | 99.03 |  | 89.14 | 91.29 |  | 51.83 | 75.76 | 82.73 | Acc
$M_{3}(0.6)$ | 0.0797 | 0.0539 |  | 0.3282 | 0.3735 |  | 1.3798 | 0.9433 | 0.5445 | Loss
98.22 | 98.31 |  | 89.25 | 87.08 |  | 53.37 | 66.68 | 82.43 | Acc
$M(0.9){\oplus}M(-0.9){\oplus}CM(0.9i)$ | 0.0873 | 0.0334 |  | 0.3671 | 0.2624 |  | 1.4029 | 0.7647 | 0.5337 | Loss
98.09 | 98.97 |  | 88.57 | 91.06 |  | 52.10 | 75.44 | 83.07 | Acc
Table 1: Loss and accuracy are reported across three different datasets, using
three distinct network architectures. The first three rows are dedicated to
benchmark optimizers, whereas the subsequent rows showcase our results. The
best benchmark result for each task (dataset and architecture pair) are
highlighted in blue. Instances where our optimizer exceeds the best baseline
result are marked in green. Additionally, the absolute best value for each
task is emphasized in bold font. The results are the average of 3 runs with
different random seeds. Figure 1:  Test accuracy graphs of RLLC and benchmark
optimizers, measured on the CIFAR-10 dataset, with the ResNet-20 network. RLLC
optimizers show faster convergence and better generalization. See related
plots and error bars in Appendix E.

RLLC based adaptive learning rate: One of the simplest case of the RLLC method
is already interesting. If there is a single memory unit containing the
momentum of previous gradients then RLLC yields an adaptive version of the
momentum SGD optimizer. In this case the learning law contains a single
coefficient, that defines an adaptively changing learning rate for the
momentum SGD. Our experiments show that this upgrade outperforms the plain
momentum SGD method, showcasing the power of RLLC. See $M(0.9)$ results in
Table 1 and on Figure 1. Note that RLLC can be applied to an arbitrary
optimizer by introducing a single memory unit storing the last learning step.
In a similar way we obtain a version of the optimizer with an adaptive
learning rate. However it may depend on the optimizer whether it leads to a
performance boost or not.

Mixing SGD and momentum SGD: We observe an intriguing phenomenon when we
enhance the memory unit of the previous method with the current gradient and
monitor the learning law throughout the training process. As shown in Figure
2, during the initial phase of training, the coefficient of the $M(0.9)$
memory unit is predominant. However, as training advances, the coefficient of
the $M(0)$ unit increases, leading to a reversal in the significance of the
two memory units. Our experiment supports Keskar and Socher (2017a) findings.
Table 1 and Figure 1 show results for $M(0.9)\oplus M(0)$. An interesting
additional detail is that in between the two extremal phases there exists a
phase which emulates the Nesterov Accelerated Gradient (NAG) method. This
occurs when the coefficient of the Momentum SGD memory unit, divided by the
coefficient of the SGD unit, equals the decay parameter of the momentum SGD.
(For more details, see Appendix A).

Figure 2:  Analysis of $M(0.9)\oplus M(0.0)$ optimizer’s memory unit’s
coefficients over time. The figure illustrates the optimizer’s transition
between momentum SGD and SGD, briefly aligning with the NAG optimizer around
the 2k step.

Multi-momentum propagators: In our experiment we investigated optimizers of
the form $M(\beta_{1})\oplus M(\beta_{2})\oplus\dots\oplus M(\beta_{k})$. We
have not optimized the hyperparameters $\beta_{i}$ but we found very promising
settings with a few trials. Our results are therefore illustrative and the
fine tuning (depending on the type of network) is subject to further research.
In Table 1 we describe our experiments with $M(0.9)\oplus M(0),M(0.9)\oplus
M(0.6),M(0.9)\oplus M(0.8)\oplus M(0.7)$ and $M(0.9)\oplus M(0.6)\oplus M(0)$
optimizers. Quite surprisingly the simplest one $M(0.9)\oplus M(0)$ (which is
mentioned earlier) is the most reliable. However on certain tasks it is
outperformed by the memory $3$ settings. Figure 3 illustrates an interesting
coupling between the coefficient of $M(0.8)$ and $M(0.7)$ memory units. See
further details in Appendix B.

Figure 3:  Analysis of the coefficients of $M(0.9)\oplus M(0.8)\oplus M(0.7)$
optimizer over time. The figure shows an interesting negative coupling between
$M(0.8)$ and $M(0.7)$. See further details in Appendix B.

$M_{m}(\lambda)$ propagator for $m\geq 2$: Jordan block propagators of the
form $M_{m}(\beta)$ and $CM_{m}(\beta)$ with $m\geq 2$ are easy to implement
in our code. In our experiments we focused on type $M_{m}(\beta)$ propagators
with $m=2,3$. It is also interesting to combine them with other propagators.
Table 1 shows results for $M_{2}(0.6)$, $M_{3}(0.6)$ and $M(0.9)\oplus
M_{2}(0.6)$. These configurations surpass the baseline optimizers in many
tasks and also surpass pure multi-momentum propagators in some specific tasks.

Complex-moment propagators: Another interesting possibility in our framework
is the usage of complex-momentum propagators. One particular example that we
experimented with is the case of $M(\beta){\oplus}M(-\beta){\oplus}CM(\beta
i)$. This choice in not random. It comes from the Jordan normal form of the
permutation matrix corresponding to the cyclic shift on $4$ elements
multiplied with $\beta$. This particular propagator is closely related to
Fourier analysis.

##  7 Limitations

Using memory units comes at a cost. Each memory unit is a vector in the
parameter space $\mathbb{R}^{n}$. In our experiments, we opted for relatively
small or medium-sized architectures. However, for architectures with a vast
parameter space, our approach with many memory units could prove to be too
memory-intensive.

It’s also worth noting that our experiments were conducted on relatively small
datasets, and future work should explore experiments on larger datasets.

##  8 Conclusion

Our experiments demonstrate that the RLLC method is capable of boosting the
performance of classical optimizers (such as SGD and momentum SGD) by
combining them and making them more adaptive. Furthermore the case of linearly
updated memory units provides a mathematically elegant framework with many new
types of promising optimizers such as the ones corresponding to larger Jordan
blocks, complex numbers and their combinations. We regard this paper as a
starting point for future research in the frame of which the full potential of
our approach is explored. One possible research direction is to introduce
adaptively changing memory update rules. In particular, in the linear setting
the pair $B\in\mathbb{R}^{k\times k},a\in\mathbb{R}^{k}$ (see Section 4) is
fixed for the whole optimization process in the current version. It would be
interesting to study a version were $B$ and $a$ are also adaptively changing
throughout learning.

##  9 Acknowledgements

The first-named author received funding from the project KPP 133921 of the
Hungarian Government. The research was also supported by the Hungarian
Ministry of Innovation and Technology NRDI Office in the framework of the
Artificial Intelligence National Laboratory Program.

## References

  * Andrychowicz et al. (2016) Marcin Andrychowicz, Misha Denil, Sergio Gomez, Matthew W. Hoffman, David Pfau, Tom Schaul, Brendan Shillingford, and Nando de Freitas.  Learning to learn by gradient descent by gradient descent, 2016.
  * Bello et al. (2017) Irwan Bello, Barret Zoph, Vijay Vasudevan, and Quoc V. Le.  Neural optimizer search with reinforcement learning, 2017.
  * Bengio et al. (1991) Y. Bengio, S. Bengio, and J. Cloutier.  Learning a synaptic learning rule.  In _IJCNN-91-Seattle International Joint Conference on Neural Networks_ , volume ii, pages 969 vol.2–, 1991.  doi: 10.1109/IJCNN.1991.155621.
  * Chen et al. (2022) Tianlong Chen, Xiaohan Chen, Wuyang Chen, Zhangyang Wang, Howard Heaton, Jialin Liu, and Wotao Yin.  Learning to optimize: a primer and a benchmark.  _J. Mach. Learn. Res._ , 23(1), jan 2022.  ISSN 1532-4435.
  * Deng (2012) Li Deng.  The mnist database of handwritten digit images for machine learning research.  _IEEE Signal Processing Magazine_ , 29(6):141–142, 2012.
  * Keskar and Socher (2017a) Nitish Keskar and Richard Socher.  Improving generalization performance by switching from adam to sgd.  12 2017a.
  * Keskar and Socher (2017b) Nitish Keskar and Richard Socher.  Improving generalization performance by switching from adam to sgd.  12 2017b.
  * Kingma and Ba (2014) Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _CoRR_ , abs/1412.6980, 2014.  URL https://api.semanticscholar.org/CorpusID:6628106.
  * Krizhevsky (2009) Alex Krizhevsky.  Learning multiple layers of features from tiny images.  2009\.  URL https://api.semanticscholar.org/CorpusID:18268744.
  * (10) Alex Krizhevsky, Vinod Nair, and Geoffrey Hinton.  Cifar-10 (canadian institute for advanced research).  URL http://www.cs.toronto.edu/~kriz/cifar.html.
  * LeCun et al. (2010) Yann LeCun, Corinna Cortes, and CJ Burges.  Mnist handwritten digit database.  _ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist_ , 2, 2010.
  * Li and Malik (2017a) Ke Li and Jitendra Malik.  Learning to optimize.  In _International Conference on Learning Representations_ , 2017a.  URL https://openreview.net/forum?id=ry4Vrt5gl.
  * Li and Malik (2017b) Ke Li and Jitendra Malik.  Learning to optimize neural nets, 2017b.
  * Li and Malik (2017c) Ke Li and Jitendra Malik.  Learning to optimize.  In _International Conference on Learning Representations_ , 2017c.  URL https://openreview.net/forum?id=ry4Vrt5gl.
  * Metz et al. (2020) Luke Metz, Niru Maheswaranathan, C. Daniel Freeman, Ben Poole, and Jascha Sohl-Dickstein.  Tasks, stability, architecture, and compute: Training more effective learned optimizers, and using them to train themselves, 2020.
  * Metz et al. (2022a) Luke Metz, C Daniel Freeman, James Harrison, Niru Maheswaranathan, and Jascha Sohl-Dickstein.  Practical tradeoffs between memory, compute, and performance in learned optimizers.  In _Conference on Lifelong Learning Agents (CoLLAs)_ , 2022a.  URL http://github.com/google/learned_optimization.
  * Metz et al. (2022b) Luke Metz, James Harrison, C. Daniel Freeman, Amil Merchant, Lucas Beyer, James Bradbury, Naman Agrawal, Ben Poole, Igor Mordatch, Adam Roberts, and Jascha Sohl-Dickstein.  Velo: Training versatile learned optimizers by scaling up, 2022b.
  * Nesterov (2012) Yurii Nesterov.  Gradient methods for minimizing composite functions.  _Mathematical Programming_ , 140:125 – 161, 2012.  URL https://api.semanticscholar.org/CorpusID:254136354.
  * Wang et al. (2022) Shipeng Wang, Yan Yang, Jian Sun, and Zongben Xu.  Variational hyperadam: A meta-learning approach to network training.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 44(8):4469–4484, 2022.  doi: 10.1109/TPAMI.2021.3061581.
  * Wichrowska et al. (2017) Olga Wichrowska, Niru Maheswaranathan, Matthew W. Hoffman, Sergio Gómez Colmenarejo, Misha Denil, Nando de Freitas, and Jascha Sohl-Dickstein.  Learned optimizers that scale and generalize.  In _Proceedings of the 34th International Conference on Machine Learning - Volume 70_ , ICML’17, page 3751–3760. JMLR.org, 2017.
  * Wilson et al. (2018) Ashia C. Wilson, Rebecca Roelofs, Mitchell Stern, Nathan Srebro, and Benjamin Recht.  The marginal value of adaptive gradient methods in machine learning, 2018\.
  * Wolpert and Macready (1997) D.H. Wolpert and W.G. Macready.  No free lunch theorems for optimization.  _IEEE Transactions on Evolutionary Computation_ , 1(1):67–82, 1997.  doi: 10.1109/4235.585893.
  * Xiao et al. (2017a) Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms, 2017a.
  * Xiao et al. (2017b) Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms, 2017b.

##  Appendix A Interpolations between SGD, momentum SGD and Nesterov
Accelerated Gradient (NAG)

As we have already explained, momentum SGD method in our interpretation is
represented as the propagator $M(\beta)$ where $0\leq\beta<1$ is the decay
parameter. In particular $M(0)$ corresponds to a memory unit which stores the
last gradient seen by the optimizer. In this sense, a memory $1$ optimizer
with memory unit $M(0)$ is basically an SGD optimizer. The learning law in
this case is a single real number which manifests as a learning rate.

If an optimizer has two memory units $M(\beta)$ and $M(0)$ then by changing
the learning law (described by a pair of real numbers) we can continuously
interpolate between momentum SGD and pure SGD. In this section we prove that
there is a third well known optimizer which can be represented in this
setting: the so-called Nesterov accelerated gradient (NAG). In the next list
we summarize the meanings of special learning laws for $M(\beta),M(0)$.

  1. 1.

$(0,c)$ : SGD

  2. 2.

$(c,0)$ : momentum SGD with decay parameter $\beta$

  3. 3.

$(c\beta,c)$ : NAG with decay parameter $\beta$

where the learning rate is $c*{\rm lr}$. Notice that since $\beta$ is a
prescribed fix number, the above three cases don’t cover the possible learning
laws for the pair $M(\beta),M(0)$.

The fact that the last learning law (if kept fixed!) is equivalent with NAG is
an easy calculation from the equations describing the NAG optimizer:

$\phi_{t+1}=\theta_{t}-r\nabla f(\theta_{t})$
$\theta_{t+1}=\phi_{t+1}+\beta(\phi_{t+1}-\phi_{t}).$

Let

$v_{t+1}:=\phi_{t}-\phi_{t+1}.$

With this notation we have the update rules:

$v_{t+1}=\beta v_{t}+r\nabla f(\theta_{t}),$ $\theta_{t+1}=\theta_{t}-(r\nabla
f(\theta_{t})+\beta v_{t+1}).$

Observe that if we introduce the update rule

$m_{t+1}:=\beta m_{t}+\nabla f(\theta_{t})$

then $v_{t}=rm_{t}$ holds at any given time $t$. Furthermore $m_{t}$
corresponds to the propagator $M(\beta)$. With this notation we have that

$\theta_{t+1}=\theta_{t}-r(\nabla f(\theta_{t})+\beta m_{t+1}).$

Assume now that $r=c*{\rm lr}$ where ${\rm lr}$ is a fixed learning rate. Then
we have that

$\theta_{t+1}=\theta_{t}-{\rm lr}(c\nabla f(\theta_{t})+c\beta m_{t+1})$

. This verifies our claim since $\nabla f(\theta_{t})$ is the propagator
$M(0)$.

##  Appendix B Additional mathematical observations

On the memory units of $M_{k}(\lambda)$ propagators: The $M_{k}(\lambda)$
propagator can naturally be interpreted as an iterated momentum propagator.
Let $m_{1},m_{2},\dots,m_{k}$ denote the the memory units. The update rule of
$M_{k}(\lambda)$ is given by

$m_{1}\leftarrow m_{1}\beta+g$ $m_{i}\leftarrow m_{i}\beta+m_{i-1}~{}{\rm
for}~{}i\geq 2.$

Thus $m_{1}$ is the momentum vector of the gradient and $m_{i}$ (for $i\geq
2$) is the momentum vector of $m_{i-1}$. One can compute that the abstract
rule of $m_{i}$ is given by the sequence
$\\{\beta^{j-i+1}{{j}\choose{i-1}}\\}_{j=0}^{\infty}.$ It follows that the
subspace generated by the abstract rules of the memory units is the space of
all sequences of the form $\\{p(j)\beta^{j}\\}_{j=0}^{\infty}$ where $p$ is a
polynomial of degree at most $k-1$. This means that we can associate such a
polynomial with each learning law. The RLLC method for $M_{k}(\lambda)$
basically adaptively navigates in this polynomial space.

Relation between Multi-momentum and $M_{k}(\lambda)$ propagators: The
progression of the learning law of $M(0.9)\oplus M(0.8)\oplus M(0.7)$ presents
an interesting phenomenon. As Figure 3 shows, the coefficients of $M(0.8)$ and
$M(0.7)$ memory units are noticeably coupled with opposite sign. One might
assume that the algorithm is just trying to cancel their effects, but the
performance improvement compared to $M(0.9)$ suggests that something more
interesting is happening here. A deeper explanation relates this optimizer to
another one of the form $M(0.9)\oplus M_{2}(0.75)$. More precisely if we
consider $\mathcal{O}(\epsilon)=M(0.9)\oplus M(0.75+\epsilon)\oplus
M(0.75-\epsilon)$ we find that as $\epsilon$ goes to $0$ the subspace spanned
by the abstract rules of the memory units converges to the subspace spanned by
the memory units of $M(0.9)\oplus M_{2}(0.75)$. In theory this convergence
means that the optimizers themselves converge. Note that in practice we can
not model $M(0.9)\oplus M_{2}(0.75)$ by $\mathcal{O}(\epsilon)$ because
numerical instability arises if $\epsilon$ is very small.

##  Appendix C Implementation details

###  C.1 Network architectures and training details

Dense network Our dense network comprises three hidden layers, each with a
width of 128 and followed by a ReLU activation function. We did not include
any normalization layers.

Convolutional network Our convolutional network features a depth of three,
with channel widths of 32, 64, and 64, each followed by a ReLU activation
function. We did not incorporate any normalization layers. Following the
convolutional layers, we apply max pooling and then a final dense layer.

ResNet-20 Our ResNet-20 variant adheres to established conventions for
CIFAR-10, employing a three-level architecture with three residual blocks at
each level. Each residual block is composed of the following sequence of
layers: Convolution-Batch Normalization-ReLU-Convolution-Batch Normalization.
A ReLU operation is applied after the addition operation in each residual
block. The convolution kernels are 3x3 in size.

###  C.2 Training details, hyperparameters

In all reported experiments, we employed a batch size of 128 and trained the
models for 10,000 iterations. We did not use a learning rate scheduler, to
avoid any potential variance in its effect across different optimizers. We run
every experiment with 3 different seed, and reported the average of the
results.

During our hyperparameter optimization process, we tested the following
potential values:

  * •

Benchmark optimizers

learning rate: 1e-7, 3e-7, 1e-6, 3e-6, 1e-5, 3e-5, 1e-4, 3e-4, 1e-3, 3e-3,
1e-2, 3e-2, 1e-1, 3e-1, 1

  * •

Our optimizers

learning rate: 0.001, 0.003, 0.01, 0.03, 0.1, 0.3

learning law - learning rate: 0.003, 0.01, 0.03

###  C.3 Trainind datasets

CIFAR-10 The CIFAR-10 dataset Krizhevsky et al.  consists of $60000$ $32x32$
colour images in $10$ classes, with $6000$ images per class. There are $50000$
training images, and $10000$ test images. We used the canonical
train–validation-test split, with $45000$ train, $5000$ validation, and
$10000$ test images. As a preprocessing, we normalized the images with the
means $(0.4914,0.4822,0.4465)$ and standard deviations
$(0.2023,0.1994,0.2010)$ for the three RGB channels, respectively. On the
ResNet task we used random resized crop (with zoom scale $0.8$-$1.2$),
horizontal flip, and random rotation.

Fashion-MNIST The Fashion-MNIST dataset Xiao et al. [2017a] consists of
$70000$ $28x28$ monochrom images in $10$ classes, with $7000$ images per
class. There are $60000$ training images, and $10000$ test images. We used the
canonical train–validation-test split, with $54000$ train, $6000$ validation,
and $10000$ test images. As a preprocessing, we normalized the images with the
mean $0.3$ and stadard deviation $0.3$.

MNIST The MNIST dataset LeCun et al. [2010] consists of $60000$ $28x28$
monocrom images in $10$ classes, with $7000$ images per class. There are
$60000$ training images, and $10000$ test images. We used the canonical
train–validation-test split, with $54000$ train, $6000$ validation, and
$10000$ test images. As a preprocessing, we normalized the images with the
mean $0.1307$ and standard deviation $0.3081$.

##  Appendix D Computational resources

For our experiments we used a server with 8 A10040GB GPUs. We reported
outcomes from a total of 16 optimizers, each optimized for hyperparameters
across seven distinct tasks (architecture - dataset pair). One round of
hyperparameter optimization with three different random seeds took
approximately 4 hours on our server, and we were able to run 16 paralelly.

Therefore, all of our results can be replicated in about 28 hours using the
same setup, or in 224 hours on a single A100 40GB GPU.

##  Appendix E Supplementary Plots

Figure 4 shows additional accuracy plots for MLP and Convolutional tasks.
Figure 5 demonstrates, that the test accuracy is consistent on different
random seeds for the demonstrated experiments.

Figure 4: Test accuracy graph of some RLLC optimizer, comparing with benchmark
optimizers. On most of the tasks RLLC optimizers perform better, than the
benchmark optimizers. Figure 5: Test accuracy graph of some RLLC optimizer,
with min-max interval, trained from $3$ random seed initialization. The
accuracy does not vary a lot, suggesting, that RLLC is robust.


# A Broadband and Compact Millimeter-Wave Imaging System based on Synthetic
Aperture Radar

Fahimeh Sepehripour1,2, Ahmad Shafiei Alavijeh2, Mohammad Fakharzadeh2,*, Amin
Khavasi2 1\. Eindhoven University of Technology, Eindhoven, Netherlands
2\. Sharif University of Technology, Tehran, Iran<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, the design, realization, and demonstration of a broadband
millimeter-wave imaging system based on the synthetic aperture radar technique
(SAR) are discussed. The proposed system, operating within the frequency range
of 25.3-30.8 GHz, consists of a tapered slot antenna as the transmitter and
two half-circle antennas as the receivers. The size of the antenna is
$19.5\times 8\,mm$ with a maximum gain of $8.5\,dB$. The transmitter and the
receiver antennas are printed on the same board. This feature leads to a
highly compact and flexible configuration, enabling the applicability of the
proposed imaging system in the handheld devices. Furthermore, it significantly
reduces the fabrication cost of the system. The proposed broadband imaging
system, being capable of performing 3D real-time imaging with high
resolutions, can be easily calibrated for each frequency within the desired
range. By performing 3D imaging from metallic objects with different shapes,
we experimentally demonstrate the high performance of the proposed system,
which offers great potentialities for a broad range of applications such as
security, medical diagnostic, concealed object detection, to name a few.

  * April 2022

##  1 Introduction

Millimeter-wave spectrum, spanning the frequency range of 30-300 GHz, has many
unique properties, making it suitable for a large variety of applications. In
particular, the electromagnetic waves within this frequency range can
penetrate through the thin layers of dielectric materials [1, 2, 3, 4, 5, 6].
This property allows one to construct the high-resolution image of many
materials such as wood, plastic and other dielectric composites that are
opaque in the visible range.

Among the different techniques proposed for millimeter-wave imaging, the so-
called synthetic aperture radar (SAR) method [7, 8, 9] is one of the most
powerful ones, allowing one to generate the image of the object in three
dimensions (3D). In this technique, a wideband single element antenna is
scanned on a two-dimensional plane located in front of the object. The signal
received by the receiver antenna is then processed using the SAR algorithms.
The major advantage of SAR is that it allows the construction of 3D images in
real-time. Apart from this, SAR enables producing images with a resolution as
high as a quarter wavelength [10, 11]. In addition, in the SAR technique, one
can achieve the image of the object by reconstructing only the real part of
the hologram.

To be able to use SAR-based imaging systems in handheld structures, it is
important to reduce their size and wight as much as possible. Mono-static
millimeter wave imaging systems represent an appropriate platform to reduce
the bulkiness of the associated structure, compared to their bi-static and
multi-static counterparts. In a mono-static imaging system, one single antenna
is used both as the transmitter and the receiver, enabling miniaturization.
Despite this advantage, mono-static imaging systems suffer from the
difficulties in the calibration process. More specifically, to properly
implement a mono-static system, it is required to calibrate the system to a
phase reference. For single frequency systems, the calibration can be readily
accomplished, since the associated phase error is the same for all points of
the synthetic aperture, whereas for a broadband system, the corresponding
phase error is not the same for all frequencies. As a result, the calibration
must be done for all frequencies individually. This renders the calibration
process extremely challenging. Nevertheless, in [12], a calibration method is
proposed, which can be conveniently performed over a wideband frequency range.
In this method, in a single board, the transmitted signal is sampled by a
receiver antenna placed at the radiating end of the transmitter antenna. The
sampled signal is then used as the reference for all frequencies, mitigating
the calibration process.

In this paper, we demonstrate a wideband millimeter wave imaging system
designed based on the SAR technique. The calibration of the proposed system is
accomplished by sampling the transmitted signal in the receiver antenna as a
reference. Since the receiver and transmitter antennas are printed on the same
board, the proposed imaging system is highly compact, enabling portability.
The antenna element of the structure is optimally designed such that the
transmitter antenna operates in a wide range of frequency with high levels of
matching. The proposed millimeter wave imaging system can have a large variety
of applications such as nondestructive testing of structures [13, 14, 15, 16],
security [17, 18, 19, 20, 21], and medical image processing [22, 23, 24].

##  2 Millimeter-Wave Imaging Modeling

###  2.1 Hardware Description

Figure 1: Conceptual sketch of the proposed millimeter wave imaging system.
The system consists of a probe antenna transmitting and receiving signals to
and from an object, a mechanical scanning system, a source, and a data
acquisitions unit.

In this section, we describe the principle of the proposed imaging system. The
goal is to realize a millimeter-wave imaging system that reconstructs the 2D
or 3D image of a target object. The proposed imaging system is based on
sending a millimeter-wave signal to the object, and analyzing the
corresponding reflected signal. Fig. 1 illustrates the proposed imaging
system. In this system, an RF signal is generated at the frequency range of
interest (25.3-30.8 GHz). The RF signal is then amplified by an amplifier
(AMMP6232) with a constant gain of 20 dB in the desired frequency range. The
corresponding amplified signal is transformed into a tapered slot antenna,
which is the transmitter antenna. The transmitter antenna sends millimeter
wave signals to the target object. It is raster-scanned in x-y direction,
where the scanning process is controlled by the computer. The corresponding
reflected field from the object is received by two receiver diodes. The based-
band signals, which are the outputs of the diodes, are then converted into two
digitial signals using a control and data acquisition unit. The associated
digital signal is then used to reconstruct the image using the generalized
synthetic aperture focusing technique (GSAFT) [25, 26, 10, 27].

###  2.2 Reconstruction Algorithm

According to the plane-wave expansion method, the received field by the
receiver antennas ($H(x^{\prime},y^{\prime})$, which is called the hologram)
is expressed as [28, 2]

$H(x^{\prime},y^{\prime})=\iint_{S}f(x,y)e^{-2jk_{0}r}dxdy,$ (1)

where $(x^{\prime},y^{\prime},0)$ and $(x,y,z_{0})$ correspond to the
locations of the antenna and the target, respectively, and $k_{0}$ is the free
space wave-number. Moreover, $f(x,y)$ corresponds to the image of the object.
In addition, $r$ is

$r=\sqrt{(x-x^{\prime})^{2}+(y-y^{\prime})^{2}+z_{0}^{2}}.$ (2)

The aim of this imaging system is to reconstruct the image of the object
$f(x,y)$. For this purpose, we decompose the exponential term in (1), which
represents a spherical wave, into a superposition of plane-waves as follows

$\begin{split}e^{-2jk_{0}r}=&e^{-2jk_{0}\sqrt{(x-x^{\prime})^{2}+(y-y^{\prime})^{2}+z_{0}^{2}}}\\\
=&\iint
e^{jk_{x^{\prime}}(x^{\prime}-x)+jk_{y^{\prime}}(y^{\prime}-y)+jk_{z}z_{0}}dk_{x^{\prime}}dk_{y^{\prime}},\end{split}$
(3)

where $k_{z}=\sqrt{4k_{0}^{2}-k_{x}^{2}-k_{y}^{2}}$. By replacing (3) in (1),
one obtains

$\begin{split}H(x^{\prime},y^{\prime})=\iint[\iint&f(x,y)e^{-jk_{x^{\prime}}x-jk_{y^{\prime}}y}dxdy]\\\
&e^{jk_{x^{\prime}}x^{\prime}+jk_{y^{\prime}}y^{\prime}+jk_{z}z_{0}}dk_{x^{\prime}}dk_{y^{\prime}},\end{split}$
(4)

in which the inner double integral represents the 2D Fourier transform of
$f(x,y)$, i.e.
$\mathcal{F}_{2D}\\{f(x,y)\\}=F(k_{x^{\prime}},k_{y^{\prime}})$. Then

$\begin{split}H(x^{\prime},y^{\prime})=\iint
F(k_{x^{\prime}},k_{y^{\prime}})e^{jk_{x^{\prime}}x^{\prime}+jk_{y^{\prime}}y^{\prime}+jk_{z}z_{0}}dk_{x^{\prime}}dk_{y^{\prime}},\end{split}$
(5)

by multiplying both sides of (5) with $e^{-jk_{z}z_{0}}$, one obtains

$\begin{split}e^{-jk_{z}z_{0}}H(x^{\prime},y^{\prime})=&\iint
F(k_{x^{\prime}},k_{y^{\prime}})e^{jk_{x^{\prime}}x^{\prime}+jk_{y^{\prime}}y^{\prime}}dk_{x^{\prime}}dk_{y^{\prime}}\\\
=&\mathcal{F}_{2D}^{-1}\\{\mathcal{F}_{2D}\\{f(x,y)\\}\\}\end{split}$ (6)

$f(x,y)$ can then be expressed according to the following equation

$f(x,y)=\mathcal{F}_{2D}^{-1}\\{\mathcal{F}_{2D}\\{H(x^{\prime},y^{\prime})e^{-jk_{z}z_{0}}\\}\\}$
(7)

$\mathcal{F}_{2D}$, and $\mathcal{F}_{2D}^{-1}$ denote the 2D Fourier
transform and inverse Fourier transform, respectively.

In the proposed system, the outputs of the diodes correspond to the real part
of $H(x^{\prime},y^{\prime})$. Here, we prove the possibility of performing
imaging using the real part of the hologram. We can rewrite the real part of
the hologram as

$\Re\\{H(x^{\prime},y^{\prime})\\}=\frac{1}{2}(H(x^{\prime},y^{\prime})+H^{*}(x^{\prime},y^{\prime}))$
(8)

in which $H(x^{\prime},y^{\prime})$ is the hologram. Based on the first term
on the right side of (8), one can obtain the image of the object at $z=z_{0}$
using the following formula [29]

$f(x,y)=\mathcal{F}_{2D}^{-1}\\{\mathcal{F}_{2D}\\{H(x^{\prime},y^{\prime})e^{-jk_{z}z_{0}}\\}\\}$
(9)

Likewise, the second term corresponds to the image of the object $z=-z_{0}$
according to the following equation

$f^{*}(x,y)=\mathcal{F}_{2D}^{-1}\\{\mathcal{F}_{2D}\\{H^{*}(x^{\prime},y^{\prime})e^{-jk_{z}(-z_{0})}\\}\\}$
(10)

This proves that, by using the real part of the hologram, we can reconstruct
the image at $z=z_{0}$.

##  3 Hardware Structure

###  3.1 Transmitter Antenna

Figure 2: Structure of the designed antenna. (a) Schematic of the top layer of the designed. (b) Bottom layer of the designed antenna. The associated dimensions are summarized in Table 1. Table 1: Dimensions of the designed Antenna Parameter | Dimension(mm)
---|---
$w_{s}$ | 8
$w_{g}$ | 0.2
$w_{m}$ | 0.5
$w_{p}$ | 0.5
$w_{r}$ | 0.2
$r_{i}$ | 0.16
$r_{o}$ | 0.31
$r_{h}$ | 1
$l_{r}$ | 2
$l_{d1}$ | 10.02
$l_{d2}$ | 5.06
$l_{1}$ | 7.47
$l_{2}$ | 12.03

A key part of a millimeter wave imaging system is the antenna element. The
antenna parameters such as the size of the antenna must be designed in a way
that the image can be reconstructed using the GSAFT algorithm, which was
described before. Fig. 2(a) and (b) show the layouts of the top and bottom
layers of the designed antenna element, respectively. The transmitter antenna
is a tapered slot antenna. The tapered slot antenna brings several advantages
[30, 31, 32]. Firstly, it enables a wideband type of operation. In addition,
it leads to a relatively high gain, low side-lobe levels, and a symmetrical
radiation pattern. The corresponding dimensions of the antenna are summarized
in Table I. The associated frequency range is 25.3-30.8 GHz, which is a
desirable region for millimeter wave imaging applications.The operation
frequency can be tuned in an arbitrary range by scaling the system dimensions.
The three-dimensional radiation pattern of the antenna is shown in Fig. 3(a),
obtained by performing direct finite-element simulations in HFSS at the
operation frequency of 26 GHz. As observed, the antenna has a symmetric
radiation pattern, in which the maximum gain ($G_{max}=8.5\,dB$) happens in
the bore-sight direction. Inset of Fig. 3(b) illustrates how the electric
field radiates from the antenna, based on which one can deduce that the
antenna has relatively wide beamwidth. Together, these features enable high-
resolution imaging based on GSAFT technique.

Figure 3: (a) Radiation pattern of the transmitter antenna. (b) Profile of the
electric field. The results are obtained at the frequency of 26 GHz.

###  3.2 Receiver

Two receiver antennas are placed symmetrically at the radiating end of the
transmitter, maintaining the radiation pattern of the antenna element
symmetric and enabling a denser spatial sampling. Each receiver is composed of
a Schottky detector diode (Skyworks SMS7621-060) and a half-circle antenna
surrounding the diode. The role of the half-circle antenna is to enhance the
coupling to the diode. The reason for choosing half-circle type antenna is its
larger electrical length compared to other elements with the same area. The
output of the Schottky diode is a real-valued DC voltage, based on which the
phase and the amplitude of the reflected electric field are extracted. This
voltage is generated by down-mixing the field scattered from the transmitter
and the associated reflected field. To characterize the performance of the
diodes, we measured the output voltage of one diode as a function of the input
power of the transmitter antenna. The associated measurement results are shown
in Fig. 4.

Figure 4: The measured output voltage of the diodes as a function of the input
power of the transmitter antenna. The measurement is performed at f=26 GHz.

###  3.3 Antenna Element Fabrication

Fig. 5(a) and (b) show the photograph of the top and bottom layers of the
fabricated antenna. The antenna is printed on a RO4003 Rogers substrate having
a thickness of 8 mils. The PCB is fabricated with a minimum gap size of 150
$\mu m$. To excite the antenna, an end-launch connector is used. The
associated S parameters are then extracted by performing scattering
measurements in an anechoic electromagnetic chamber. Fig. 5(c) represents the
measured (dashed blue line) scattering parameter $(S_{11})$ and compares it
with simulation (solid red line). The simulation result is obtained by
performing the full-wave simulation in Ansys HFSS. As observed, within the
frequency range of 25.3-30.8 GHz, $S_{11}$ of the antenna is below 10 dB
(purple line).

Figure 5: Performance analysis of the fabricated antenna. (a), and (b)
Photographs of the top (panel a) and bottom (panel b) of the fabricated
antenna. (c) The corresponding measured (blue line) and simulated (red line)
$S_{11}$.

###  3.4 Data Acquisition System

To analyze the information received by the antennas, we designed a controller
shown in Fig. 6(a). The controller includes an amplifier, an offset-removal
circuit and an analog to digital converter (ADC). Fig. 6(b) depicts the
schematic of the offset-removal circuit. As observed, it includes a non-
inverting operational amplifier (LM224D) and a potentiometer (left panel),
used to maximize the amplification. The dynamic range of the output is
amplified using an operational amplifier, as long as it is not saturated. The
values of the circuit elements are shown in Fig. 6(b). For realizing the ADC,
we have used an ATmega328 micro-controller in Arduino-nano. The digital signal
is transferred to a computer using USB connectors. The corresponding digital
signal is used to reconstruct the image using the GSAFT algorithm described in
section II.

Figure 6: Photograph of the data acquisition system. (a) the controller
consists of Arduino-nano, LM224D amplifiers, and an offset-removal circuit.
(b) Schematic of the off-set removal circuit of our designed controller

##  4 Experimental Results

To verify the performance of the proposed imaging system, we used a metallic
gun model as the target object. Fig. 7(a) shows the experimental setup, where
the antenna element is raster-scanned over a rectangular plane (note that the
center of this rectangular plane is matched with that of the object. The
distance between the center of antenna and that of the object is 39 mm). For
mechanical scanning, the step distance along horizontal and vertical
directions are 2.5 mm and 0.5 mm, respectively. The imaging is performed at a
single frequency ($f=26\,GHz$). Fig. 7(b) and (c) show the corresponding
holograms based on the output voltages of the upper and lower diodes,
respectively. As it is observed, in both cases, the image of the object is
properly reconstructed.

Figure 7: Experimental measurements. (a) Photograph of the experimental setup
used to perform millimeter wave imaging. (b), (c) Corresponding images of the
gun, based on the signal received at the antenna 1 and 2. Figure 8: Millimeter
wave imaging of a metallic cross-shaped object. (a) Photograph of a metallic
cross. (b) Reconstructed image of the object. (c), (d) Corresponding hologram
associated with the antenna 1 at the frequency of 26.1 GHz and 28.1 GHz,
respectively.

In order to perform wideband imaging, we fed the antenna element by a wideband
millimeter-wave source in the frequency range of 25.1-30.8 GHz (the frequency
step is 0.5 GHz). We then performed imaging of a cross-shaped object, whose
photograph is shown in Fig. 8(a). Fig. 8(b) shows the corresponding
reconstructed image. As observed, the wideband imaging system allows us to
reconstruct the sharp edges of the image more precisely. Furthermore, wideband
imaging enables focusing on depth locations. This feature is clear from the
result of Fig. 8(b), where the reflected field has almost the same amplitude
inside the cross area. Insets of Fig. 8(c) ,(d) represent the corresponding
hologram associated with the first antenna at two different frequencies,
namely 26.1 GHz (panel c) and 28.1 GHz (panel d).

##  5 Conclusion

We designed and experimentally demonstrated a wideband millimeter wave imaging
system based on SAR technique. By performing measurements, we showed that the
proposed imaging system was capable of reconstructing three-dimensional images
of arbitrary shaped objects. Compared to traditional microwave imaging systems
in which the transmitter and receiver antennas are often printed on two
different boards, both transmitter and receiver antennas of our proposed
system were printed on the same board. This provides us with a unique
opportunity to tackle the problem of phase referencing for different
frequencies. The proposed imaging system has potential applications in
security scanning, non-destructive sensing, spectroscopy, and medical
screening. By substituting the single element antenna of our proposed
structure with a linear array, one can benefit from the advantages of
electrical scanning.

## Reference

## References

  * [1] J. Moll, P. Schops, and V. Krozer, “Towards three-dimensional millimeter-wave radar with the bistatic fast-factorized back-projection algorithm—potential and limitations,” _IEEE Transactions on Terahertz Science and Technology_ , vol. 2, no. 4, pp. 432–440, 2012.
  * [2] D. M. Sheen, D. L. McMakin, and T. E. Hall, “Three-dimensional millimeter-wave imaging for concealed weapon detection,” _IEEE Transactions on microwave theory and techniques_ , vol. 49, no. 9, pp. 1581–1592, 2001.
  * [3] A. H. Lettington, M. R. Yallop, and D. Dunn, “Review of super-resolution techniques for passive millimeter-wave imaging,” in _Infrared and Passive Millimeter-wave Imaging Systems: Design, Analysis, Modeling, and Testing_ , vol. 4719. International Society for Optics and Photonics, 2002, pp. 230–239.
  * [4] C. Jia, G. Geng, T. Chang, C. Yang, W. Fan, Y. Guo, Z. Sun, L. Liu, Q. Guo, and H.-L. Cui, “Millimeter wave imaging of cracks in bulk coal,” _IEEE Transactions on Terahertz Science and Technology_ , vol. 6, no. 4, pp. 554–562, 2016.
  * [5] V. M. Patel, J. N. Mait, D. W. Prather, and A. S. Hedden, “Computational millimeter wave imaging: problems, progress, and prospects,” _IEEE Signal Processing Magazine_ , vol. 33, no. 5, pp. 109–118, 2016.
  * [6] R. Appleby, D. A. Robertson, and D. Wikner, “Millimeter wave imaging: a historical review,” in _Passive and Active Millimeter-Wave Imaging XX_ , vol. 10189. International Society for Optics and Photonics, 2017, p. 1018902.
  * [7] A. Elboushi and A. Sebak, “B14. active millimeter-wave imaging system for hidden weapons detection,” in _2012 29th National Radio Science Conference (NRSC)_. IEEE, 2012, pp. 111–118.
  * [8] S. Camlica, A. C. Gurbuz, and O. Arikan, “Autofocused spotlight sar image reconstruction of off-grid sparse scenes,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 53, no. 4, pp. 1880–1892, 2017.
  * [9] H. Zamani, M. Fakharzadeh, A. Amini, and F. Marvasti, “Random sampling with iterative recovery for millimeter-wave imaging systems,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , 2020.
  * [10] K. Qin, C. Yang, and F. Sun, “Generalized frequency-domain synthetic aperture focusing technique for ultrasonic imaging of irregularly layered objects,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 61, no. 1, pp. 133–146, 2014.
  * [11] M. T. Ghasr, M. J. Horst, M. R. Dvorsky, and R. Zoughi, “Wideband microwave camera for real-time 3-d imaging,” _IEEE Transactions on Antennas and Propagation_ , vol. 65, no. 1, pp. 258–268, 2016.
  * [12] M. T. Ghasr, J. T. Case, and R. Zoughi, “Novel reflectometer for millimeter-wave 3-d holographic imaging,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 63, no. 5, pp. 1328–1336, 2014.
  * [13] M. Kempin, M. T. Ghasr, J. T. Case, and R. Zoughi, “Modified waveguide flange for evaluation of stratified composites,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 63, no. 6, pp. 1524–1534, 2013.
  * [14] S. Agarwal and D. Singh, “An adaptive statistical approach for non-destructive underline crack detection of ceramic tiles using millimeter wave imaging radar for industrial application,” _IEEE Sensors Journal_ , vol. 15, no. 12, pp. 7036–7044, 2015.
  * [15] K. M. Shafi, M. S. ur Rahman, M. A. Abou-Khousa, and M. R. Ramzi, “Applied microwave imaging of composite structures using open-ended circular waveguide,” in _2017 IEEE International Conference on Imaging Systems and Techniques (IST)_. IEEE, 2017, pp. 1–5.
  * [16] H. Alidoustaghdam and M. Çayören, “Non-destructive testing of concrete tunnels with qualitative microwave imaging,” in _2020 German Microwave Conference (GeMiC)_. IEEE, 2020, pp. 29–31.
  * [17] R. Appleby and R. N. Anderton, “Millimeter-wave and submillimeter-wave imaging for security and surveillance,” _Proceedings of the IEEE_ , vol. 95, no. 8, pp. 1683–1690, 2007.
  * [18] D. M. Sheen, D. L. McMakin, and T. E. Hall, “Active millimeter-wave and sub-millimeter-wave imaging for security applications,” in _2011 International Conference on Infrared, Millimeter, and Terahertz Waves_. IEEE, 2011, pp. 1–3.
  * [19] Z. Wu, L. Zhang, H. Liu, and N. Kou, “Enhancing microwave metamaterial aperture radar imaging performance with rotation synthesis,” _IEEE Sensors Journal_ , vol. 16, no. 22, pp. 8035–8043, 2016.
  * [20] R. Zhang and S. Cao, “Real-time human motion behavior detection via cnn using mmwave radar,” _IEEE Sensors Letters_ , vol. 3, no. 2, pp. 1–4, 2018.
  * [21] Z. Briqech, S. Gupta, A. Beltayib, A. Elboushi, A.-R. Sebak, and T. A. Denidni, “57-64 ghz imaging/detection sensor–part i: System setup and experimental evaluations,” _IEEE Sensors Journal_ , vol. 20, no. 18, pp. 10 824–10 832, 2020.
  * [22] T. Henriksson, N. Joachimowicz, C. Conessa, and J.-C. Bolomey, “Quantitative microwave imaging for breast cancer detection using a planar 2.45 ghz system,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 59, no. 10, pp. 2691–2699, 2010.
  * [23] M. Klemm, J. A. Leendertz, D. Gibbins, I. J. Craddock, A. Preece, and R. Benjamin, “Microwave radar-based differential breast cancer imaging: Imaging in homogeneous breast phantoms and low contrast scenarios,” _IEEE Transactions on Antennas and Propagation_ , vol. 58, no. 7, pp. 2337–2344, 2010.
  * [24] M. A. Aldhaeebi, T. S. Almoneef, H. Attia, and O. M. Ramahi, “Near-field microwave loop array sensor for breast tumor detection,” _IEEE Sensors Journal_ , vol. 19, no. 24, pp. 11 867–11 872, 2019.
  * [25] M. Rezaei, H. Zamani, M. Fakharzadeh, and M. Memarian, “Quality improvement of millimeter-wave imaging systems using optimized dual polarized arrays,” _IEEE Transactions on Antennas and Propagation_ , vol. 69, no. 10, pp. 6848–6856, 2021.
  * [26] H. Zamani and M. Fakharzadeh, “1.5-d sparse array for millimeter-wave imaging based on compressive sensing techniques,” _IEEE Transactions on Antennas and Propagation_ , vol. 66, no. 4, pp. 2008–2015, 2018.
  * [27] E. Nilsson and L. Baath, “Radar interferometric measurements with a planar patch antenna array,” _IEEE Sensors Journal_ , vol. 7, no. 7, pp. 1025–1031, 2007.
  * [28] M. Soumekh, “Bistatic synthetic aperture radar inversion with application in dynamic object imaging,” in _Acoustics, Speech, and Signal Processing, IEEE International Conference on_. IEEE Computer Society, 1991, pp. 2577–2578.
  * [29] D. M. Sheen, “Sparse multi-static arrays for near-field millimeter-wave imaging,” in _2013 IEEE Global Conference on Signal and Information Processing_. IEEE, 2013, pp. 699–702.
  * [30] K. Ebnabbasi, D. Busuioc, R. Birken, and M. Wang, “Taper design of vivaldi and co-planar tapered slot antenna (tsa) by chebyshev transformer,” _IEEE transactions on antennas and propagation_ , vol. 60, no. 5, pp. 2252–2259, 2012\.
  * [31] M. A. Abou-Khousa, M. T. Ghasr, S. Kharkovsky, D. Pommerenke, and R. Zoughi, “Modulated elliptical slot antenna for electric field mapping and microwave imaging,” _IEEE transactions on antennas and propagation_ , vol. 59, no. 3, pp. 733–741, 2010.
  * [32] S. Sugawara, Y. Maita, K. Adachi, K. Mori, and K. Mizuno, “A mm-wave tapered slot antenna with improved radiation pattern,” in _1997 IEEE MTT-S International Microwave Symposium Digest_ , vol. 2. IEEE, 1997, pp. 959–962.