# Wave Physics-Informed Matrix Factorizations

Harsha Vardhan Tetali, Joel B. Harley, Benjamin D. Haeffele

###### Abstract

With the recent success of representation learning methods, which includes
deep learning as a special case, there has been considerable interest in
developing techniques that incorporate known physical constraints into the
learned representation. As one example, in many applications that involve a
signal propagating through physical media (e.g., optics, acoustics, fluid
dynamics, etc), it is known that the dynamics of the signal must satisfy
constraints imposed by the wave equation. Here we propose a matrix
factorization technique that decomposes such signals into a sum of components,
where each component is regularized to ensure that it nearly satisfies wave
equation constraints. Although our proposed formulation is non-convex, we
prove that our model can be efficiently solved to global optimality. Through
this line of work we establish theoretical connections between wave-informed
learning and filtering theory in signal processing. We further demonstrate the
application of this work on modal analysis problems commonly arising in
structural diagnostics and prognostics.

###### Index Terms:

physics-guided machine learning, dictionary learning, modal decomposition,
unsupervised learning

††publicationid: pubid:  © 2023 IEEE, DOI: 10.1109/TSP.2023.3348948

##  I Introduction

Representation learning has gained importance in fields that utilize data
generated through physical processes, such as weather forecasting [1],
manufacturing [2], structural health monitoring [3], acoustics [4], and
medical imaging [5]. Moreover, learning physically consistent and
interpretable features can improve our understanding of the physically viable
information about the data and the composition of the system or process
generating it. Unfortunately, however, features of the data learned through
generic machine learning algorithms typically do not correspond to physically
interpretable quantities.

For example, physically consistent solutions are important to spatio-temporal
modal analysis, where recovering modes of vibration yield important
information about a physical system or process (e.g., buildings, bridges,
space shuttles, satellites, acoustic instruments) [6, 7, 8, 9, 10].
Furthermore, these problems are also relevant to frequency-wavenumber spectral
analysis [11], which is the basis for many advances in RADAR [11], seismic
[12], and acoustic [13] beamforming and signal processing, where each sinusoid
may correspond to a different source or target. Physically consistent
solutions provide a trustworthy basis in terms of the physics of the system in
contrast to any other basis, which might be obtained by generic machine
learning algorithms. This suggests the need for learning algorithms which
choose bases that are physically consistent.

###  I-A Related Work

Modified learning paradigms, suited to different physical application domains,
have begun to draw interest [14]. For example, several researchers have
designed physics-informed neural networks to learn approximate solutions to a
partial differential equation [15, 16, 17, 18, 19]. These have been recently
applied to ultrasonic surface waves to extract velocity parameters and denoise
data [20]. Similar physics-guided neural networks [14, 21, 22, 23] use
physics-based regularization to improve data processing, demonstrating that
regression tasks (such as estimating the temperature throughout a lake) can be
more accurate and robust when compared with purely physics-based solutions or
purely data-driven solutions. These approaches demonstrate the strong
potential for physics-informed machine learning but remain early in their
study.

Representation theory has been applied to other wave-based problems. For
example, generative models have been studied for applications in
electroencephalography [24] and seismology [25], although these learning
systems do not assume physical knowledge. Other prior work has considered
sparse generative representation of ultrasonic waves through the use of sparse
signal processing [26], dictionary learning [27, 28, 29], and neural network
sparse autoencoders [30]. These methods utilize the fact that many waves can
be represented as a small, sparse sum of spatial modes and attempt to extract
these modes within the data [27]. These dictionary learning algorithms have
further been combined with wave-informed regularizers to incorporate physical
knowledge of wave propagation into the solution [31, 32].

One of the key challenges for many of these problems is that they typically
require solving a non-convex optimization problem, and as a result the learned
representation can be highly dependent on details such as initialization or
choice of optimization algorithm. Attention has been devoted recently to the
non-convex optimization problems that arise in representation learning, and
positive results have begun to emerge in certain settings. For example, within
the context of low-rank matrix recovery, it has been shown that gradient
descent is guaranteed to converge to a local minimum and all local minima will
be globally optimal under certain circumstances [33, 34, 35, 36]. Further,
other recent work has studied ‘structured’ matrix factorization problems,
where one promotes properties in the matrix factors by imposing regularization
on the factorized matrices (e.g., if one desires a factor to be sparse an
$\ell_{1}$ norm penalty could be imposed on one of the matrix factors). It has
been shown that such problems can also be solved to global optimality in
certain circumstances, but whether such guarantees can be made depends
critically on how the regularization is formulated [37, 38, 39].

###  I-B Contributions

In this work we make the following contributions. We _first_ introduce a wave-
informed matrix factorization model, which ensures that the recovered
components minimize wave equation constraints, allowing for direct
interpretability of the factors based on assumed physics. This is
distinguished from traditional factorization methods, such as principal
component analysis (PCA), which are not constrained by physical knowledge as
well as from models which explicitly parameterize wave-like factors, and are
often too constrained for a problem, particularly in the presence of noise. In
other words, our approach balances the utilization of both prior assumptions
from physics and fills the gap between purely data driven methods [40] and
purely physics based methods.

_Second_ , although our resulting formulation requires solving a _non-convex_
optimization problem, we identify its form amenable to be provably solved to
global optimality [38]. This differs from prior wave-physics informed
dictionary learning [31], which provides no such guarantees.

_Third_ , unlike physics-informed neural networks, we show that our wave
physics-informed matrix factorization has a satisfying interpretation as an
optimal filter problem. Specifically, we show that our optimization algorithm
based on [38] reveals a solution is that of an optimal bank of 1st-order
Butterworth filters. This interpretability enables confidence in the algorithm
and intelligent selection of hyper-parameters.

In summary, this work’s novel contributions are in demonstrating the first
framework for matrix factorization that is constrained solely on the wave
equation. Unlike wave-informed dictionary learning [31], we show this
algorithm is both highly interpretable from a signal processing perspective
and satisfies global optimality conditions based on results the literature
[38]. Finally, we evaluate the effectiveness of the algorithm in four
scenarios – homogeneous vibrations, non-homogeneous vibrations, travelling
plane waves, and segmented vibrations. We show that this framework outperforms
most other methods used the literature for wave mode decomposition,
particularly in the presence of high noise and segmented, spatially varying
wave behavior.

##  II Wave-Informed Matrix Factorization

In this section, we develop a novel matrix factorization framework that allows
us to decompose a matrix $\mathbf{Y}$, with each column representing sampling
along one dimension (e.g., space), as the product of two matrices
$\mathbf{D}\mathbf{X}^{\top}$, where one of the matrices will be softly
constrained to nearly satisfy wave equation constraints. Additionally, our
model will learn the number of modes (columns in $\mathbf{D}$ and
$\mathbf{X}$) directly from the data without specifying this a priori.

In subsection II-A, we introduce our objective function. This subsection can
be seen as an extension of the previously proposed wave-informed dictionary
learning framework to a matrix factorization framework. In II-B, we go over a
few ideas from the optimization literature that are useful to solve the
formulated problem. In II-C and II-D, we describe the main portions of the
wave-physics informed matrix factorization algorithm. In II-E, we describe an
interesting connection of this framework to signal processing. We emphasize
that the novel parts of this work is in subsections II-C and II-E. The work
optimization discussions in II-B are specific cases of previous work and are
shown for the benefit of the reader, particularly those in wave physics who
may not be familiar with the optimization literature.

###  II-A Problem Formulation

The problem formulation described in this section is inspired from the wave-
physics regularizer proposed in [31]. In this work, we additionally justify
the wave-physics informed problem formulation in [31] as a naturally arising
idea from the theory of solutions of partial differential equations. A key
concept in the solution of partial differential equations is the notion of
separation. For example, given a PDE in space $\ell$ and time $t$, we can
often assume the solution is of the form

$\displaystyle f(\ell,t)=d(\ell)x(t)\;.$ (1)

In the case of linear PDEs, we find that if $d_{i}(\ell)x_{i}(t)$ are
solutions for all $i\in[1,N]$, then

$\displaystyle f_{i}(\ell,t)=\sum_{i=1}^{N}{\alpha_{i}}d_{i}(\ell)x_{i}(t)$
(2)

is also a solution for any arbitrary choice of linear coefficients,
$\\{\alpha_{i}\\}_{1=1}^{N}$. For a discrete approximation of this solution,
we note that the product of continuous functions $d_{i}(\ell)x_{i}(t)$ on
$\ell\in[0,L]$ and $t\in[0,T]$ can be approximated as an outer product of two
vectors, leading directly to a matrix factorization model, where we further
note that without loss of generality the scaling coefficients $(\alpha_{i})$
can be absorbed into the factorized matrices.

Specifically, we consider matrices $(\mathbf{D}\in\mathbb{R}^{N_{d}\times
N},\mathbf{X}\in\mathbb{R}^{N_{t}\times N})$ defined such that

$\displaystyle\mathbf{d}_{i}$
$\displaystyle=\left[d_{i}(0),d_{i}(\Delta\ell),\cdots,d_{i}(N_{d}\Delta\ell)\right]^{\top}$
$\displaystyle\mathbf{x}_{i}$ $\displaystyle=\left[x_{i}(0),x_{i}(\Delta
t),\cdots,x_{i}(N_{t}\Delta t)\right]^{\top}$
$\displaystyle\mathbf{D}\mathbf{X}^{\top}$
$\displaystyle=\sum_{i=1}^{N}d_{i}(\ell)x_{i}(t)=\sum_{i=1}^{N}\mathbf{d}_{i}\mathbf{x}_{i}^{\top}\;,$
(3)

where $\mathbf{d}_{i}$ and $\mathbf{x}_{i}$ denote columns of $\mathbf{D}$ and
$\mathbf{X}$, respectively, and
$N_{d}=\left\lfloor{L/\Delta\ell}\right\rfloor$ and
$N_{t}=\left\lfloor{T/\Delta t}\right\rfloor)$. Recall, we also wish to
enforce physical consistency in our model, which we accomplish through theory-
guided regularization (see [14]) in a cost function which 1) captures how our
matrix factorization model matches the observed samples and 2) incorporates a
regularization term which enforces physical consistency in the matrix factors.
Specifically, we consider a model of the form

$\underset{\mathbf{D},\mathbf{X}}{\text{ }\min\text{
}}{\tfrac{1}{2}\|\mathbf{Y}-\mathbf{D}\mathbf{X}^{\top}\|^{2}_{F}+\lambda\Theta(\mathbf{D},\mathbf{X})}$
(4)

where $\|\mathbf{Y}-\mathbf{D}\mathbf{X}^{\top}\|^{2}_{F}$ is a data
consistency loss, $\Theta(\mathbf{D},\mathbf{X})$ is our physics-informed
regularizer, and $\lambda>0$ is a hyper-parameter to balance the trade-off
between fitting the data and satisfying our model assumptions. When
$\Theta(\mathbf{D},\mathbf{X})$ is defined by the wave equation, we call this
_wave-informed matrix factorization_.

We derive a wave-informed $\Theta(\mathbf{D},\mathbf{X})$ by considering the
one-dimensional wave equation evaluated at $d_{i}(l)x_{i}(t)$

$\displaystyle\frac{\partial^{2}\left[d_{i}(\ell)x_{i}(t)\right]}{\partial
l^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{c^{2}}\frac{\partial^{2}\left[d_{i}(\ell)x_{i}(t)\right]}{\partial
t^{2}}$ (5)

The above equation constrains that each $d_{i}(\ell)x_{i}(t)$ component
($i\in[1,N]$) is a wave. This can be seen as decomposing the data into a
superposition of simpler waves. Furthermore, the Fourier transform in time on
both sides of (5) yields

$\displaystyle\frac{\partial^{2}\left[d_{i}(\ell)\right]}{\partial
l^{2}}X_{i}(\omega)$ $\displaystyle=$
$\displaystyle\frac{-\omega^{2}}{c^{2}}d_{i}(\ell)X_{i}(\omega)$ (6)
$\displaystyle\frac{\partial^{2}\left[d_{i}(\ell)\right]}{\partial\ell^{2}}$
$\displaystyle=$ $\displaystyle\frac{-\omega^{2}}{c^{2}}d_{i}(\ell)$ (7)

at points where $X_{i}(\omega)\neq 0$. This is specifically known as the
Helmholtz equation. The Helmholtz equation is analytically solvable for
constant $\omega/c$ and the solution takes the form
$d_{i}(\ell)=A\sin(\omega\ell/c)+B\cos(\omega\ell/c)$, and the constants $A$
and $B$ can be determined by initial conditions. While the Helmholtz equation
is often applied in the frequency domain, it is also applicable in time, where
the relationship must be true for every frequency-velocity ($\omega,c$) mode
in the data. Hence, if the data is sparse in these domains (as is often the
case in modal analysis), then we can use the Helmholtz constraint to estimate
and decompose all ($\omega,c$) modes in $\mathbf{Y}$.

We discretize the Helmholtz equation such that

$\displaystyle\mathbf{L}\mathbf{d}_{i}$ $\displaystyle=$ $\displaystyle-
k_{i}^{2}\mathbf{d}_{i}$ (8)

where $k_{i}=\omega_{i}/c_{i}$, also known as the angular wavenumber, and
$\mathbf{L}$ is a discrete second derivative operator. The specific form of
$\mathbf{L}$ depends on the boundary conditions (e.g., Dirichlet, Dirichlet-
Neumann, etc), which can be readily found in the numerical methods literature
(see for example [41, 42]). In this paper, we use the Laplacian matrix
($\mathbf{L}$) defined by

$\displaystyle\mathbf{L}$ $\displaystyle=$ $\displaystyle\frac{1}{(\Delta
l)^{2}}\begin{bmatrix}-2&1&0&0&0&\cdots&0\\\ 1&-2&1&0&0&\cdots&0\\\
0&1&-2&1&0&\cdots&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&0&0&\cdots&-2\\\ \end{bmatrix}.$ (9)

Note that to reduce notational complexity, we assume $\Delta l=1$ without any
loss of generality, as any change in $\Delta l$ can be absorbed into the value
of $k_{i}$ after plugging (9) into (8).

We also specify the dependence on $i$ owing to the fact that $\omega/c$ and
$d_{i}(l)$ determine each other. For the factorization to satisfy the wave
equation, we desire that (8) is satisfied for some value of $k_{i}$, which we
promote with the regularizer

$\displaystyle\min_{k_{i}}\|\mathbf{L}\mathbf{d}_{i}+k_{i}^{2}\mathbf{d}_{i}\|^{2}_{F}\;.$
(10)

Note that for the above regularizer it is easily shown that the optimal value
for $k_{i}^{2}$ must lie in the range $[0,\frac{4}{\Delta l}]$, where
$\frac{4}{\Delta l}$ is the smallest eigenvalue of $-\mathbf{L}$ [43]. Since
we specifically chose $\Delta l=1$ (without loss of generality), we have that
the optimal value of $k_{i}^{2}$ must lie in the range $[0,4]$.

In addition to requiring that the columns of $\mathbf{D}$ satisfy the wave
equation, we also would like to constrain the number of modes to be minimal.
We accomplish this by adding the squared Frobenius norms to both $\mathbf{D}$
and $\mathbf{X}$, which is known to induce low-rank solutions in the product
$\mathbf{D}\mathbf{X}^{\top}$ due to connections with the variational form of
the nuclear norm [38, 37, 44]. Hence, our complete regularization is defined
by

$\displaystyle\Theta(\mathbf{D},\mathbf{X},\mathbf{k})$
$\displaystyle=\tfrac{1}{2}\sum_{i=1}^{N}\bar{\theta}(\mathbf{d}_{i},\mathbf{x}_{i},k_{i})$
(11) $\displaystyle\bar{\theta}(\mathbf{d}_{i},\mathbf{x}_{i},k_{i})$
$\displaystyle=\left(\|\mathbf{x}_{i}\|_{F}^{2}+\|\mathbf{d}_{i}\|_{F}^{2}\right)+\\!\gamma\min_{k_{i}}\|\mathbf{L}\mathbf{d}_{i}\\!+\\!k_{i}^{2}\mathbf{d}_{i}\|_{F}^{2}$
$\displaystyle=\|\mathbf{x}_{i}\|_{F}^{2}+\min_{k_{i}}\mathbf{d}_{i}^{\top}\left(\mathbf{I}+\gamma(\mathbf{L}+\mathbf{I}k_{i}^{2})^{2}\right)\mathbf{d}_{i}$
$\displaystyle=\|\mathbf{x}_{i}\|_{F}^{2}+\min_{k_{i}}\mathbf{d}_{i}^{\top}\mathbf{A}(k_{i})\mathbf{d}_{i}$
(12) $\displaystyle\text{where}\ \mathbf{A}(k_{i})$
$\displaystyle=\mathbf{I}+\gamma(\mathbf{L}+\mathbf{I}k_{i}^{2})^{2}.$ (13)

With this regularization function, the complete cost function is defined by

$\displaystyle\min_{\begin{subarray}{c}\mathbf{D},\mathbf{X},\mathbf{k},N\end{subarray}}\\!\tfrac{1}{2}\|\mathbf{Y}-\mathbf{D}\mathbf{X}^{\top}\|_{F}^{2}\\!+\\!\tfrac{\lambda}{2}\\!\|\mathbf{X}\|_{F}^{2}\\!+\\!\tfrac{\lambda}{2}\\!\sum_{i=1}^{N}\mathbf{d}_{i}^{\top}\mathbf{A}(k_{i})\mathbf{d}_{i}$
(14)

where note that we are also optimizing over the number of modes $(N)$ and the
wavenumber $\mathbf{k}=\left[k_{1},\cdots,k_{N}\right]^{\top}$ for each mode
in the data.

###  II-B Model Optimization

Our model in (14) is inherently non-convex in
$(\mathbf{D},\mathbf{X},\mathbf{k})$ due to the matrix factorization model and
further complicated by the fact that we are additionally searching over the
number of columns/modes $N$. However, despite the challenge of non-convex
optimization, we can solve (14) by leveraging prior results from optimization
theory for structured matrix factorization [37, 39]. Note that, for the
remainder of this subsection we reiterate the concepts and framework described
in [37, 38] leading to adapting the algorithm mentioned therein to our
specific model and finally obtaining an optimal solution of the objective
function (14). For a matrix factorization problem of the form

$\min_{N}{\min_{\mathbf{D},\mathbf{X}}}\;L(\mathbf{D}\mathbf{X}^{\top})+\lambda\sum_{i=1}^{N}\bar{\theta}(\mathbf{d}_{i},\mathbf{x}_{i})$
(15)

where $L(\widehat{\mathbf{Y}})$ is any function which is convex and once
differentiable in $\widehat{\mathbf{Y}}$ and
$\widebar{\theta}(\mathbf{d},\mathbf{x})$ is any function that satisfies the
following three conditions

  1. 1.

$\widebar{\theta}(\alpha\mathbf{d},\alpha\mathbf{x})=\alpha^{2}\widebar{\theta}(\mathbf{d},\mathbf{x}),\
\forall(\mathbf{d},\mathbf{x})$, $\forall\alpha\geq 0$.

  2. 2.

$\widebar{\theta}(\mathbf{d},\mathbf{x})\geq 0,\
\forall(\mathbf{d},\mathbf{x})$.

  3. 3.

For all sequences $(\mathbf{d}^{(n)},\mathbf{x}^{(n)})$ such that
$\|\mathbf{d}^{(n)}(\mathbf{x}^{(n)})^{\top}\|\rightarrow\infty$ then
$\widebar{\theta}(\mathbf{d}^{(n)},\mathbf{x}^{(n)})\rightarrow\infty$ ,

an algorithm exists that can obtain the global minimum solution. In
supplementary material VII-A, we show that the wave-informed cost function
satisfies all three conditions, where the optimization over the $k_{i}$
parameters is done inside $\theta$ as in (II-A).

Moreover, in [38] it is shown that a given point
$(\widetilde{\mathbf{D}},\widetilde{\mathbf{X}})$ is a globally optimal
solution of (15) iff the following two conditions are satisfied:

  1. 1.

$\langle-\nabla
L(\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top}),\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top}\rangle=\lambda\sum_{i=1}^{N}\bar{\theta}(\widetilde{\mathbf{d}}_{i},\widetilde{\mathbf{x}}_{i})$

  2. 2.

$\Omega_{\bar{\theta}}^{\circ}(-\tfrac{1}{\lambda}\nabla
L(\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top}))\leq 1$

where in our scenario

$\displaystyle-\tfrac{1}{\lambda}\nabla(L(\mathbf{D}\mathbf{X}^{\top}))$
$\displaystyle=\tfrac{1}{\lambda}\left(\mathbf{Y}-\mathbf{D}\mathbf{X}^{\top}\right)$

denotes the gradient w.r.t. the matrix product and
$\Omega_{\bar{\theta}}^{\circ}(\cdot)$ is referred to as the polar problem,
which is defined as

$\Omega_{\bar{\theta}}^{\circ}(\mathbf{Z})\equiv\sup_{\mathbf{d},\mathbf{x}}\mathbf{d}^{\top}\mathbf{Z}\mathbf{x}\
\ \textnormal{s.t.}\ \ \bar{\theta}(\mathbf{d},\mathbf{x})\leq 1.$ (16)

It is further shown in [38] (Proposition 3) that the first condition above
will always be satisfied for any first-order stationary point
$(\widetilde{\mathbf{D}},\widetilde{\mathbf{X}})$, and that if a given point
is not globally optimal then the objective function (15) can always be
decreased by augmenting the current factorization by a solution to the polar
problem as a new column:

$\displaystyle(\mathbf{D},\mathbf{X})$
$\displaystyle\leftarrow\left(\left[\widetilde{\mathbf{D}},\
\tau\mathbf{d}^{*}\right],\left[\widetilde{\mathbf{D}},\
\tau\mathbf{x}^{*}\right]\right)\\!$ (17)
$\displaystyle\\!\mathbf{d}^{*},\mathbf{x}^{*}$
$\displaystyle\in\operatorname*{arg\,max}_{\mathbf{d},\mathbf{x}}\tfrac{1}{\lambda}\mathbf{d}^{\top}\\!\left(\mathbf{Y}-\mathbf{D}\mathbf{X}^{\top}\right)\mathbf{x}\
\ \textnormal{s.t.}\ \ \bar{\theta}(\mathbf{d},\mathbf{x})\leq 1$

for an appropriate choice of step size $\tau>0$.

Despite the fact that [38] gives a meta-algorithm that is guaranteed to solve
(15) to global optimality, one needs to solve the polar problem in (16) to
both verify global optimality and escape non-optimal stationary points, which
is itself another optimization problem. This presents a potentially
significant computational challenge in practice, as for many choices of
regularization function, solving the polar problem is NP-hard [39]. For
example, if the regularization on the matrix columns is simply the sum of the
$\ell_{2}$ norms of the columns, then the polar problem is easily solved as a
maximum eigenvalue problem, but if we either replace one of the $\ell_{2}$
norms by the $\ell_{\infty}$ norm, the problem becomes NP-hard [39].
Fortunately, however, in the next section we show that the polar problem that
results from our wave-informed factorization model is tractable and leads to
polynomial time global optimality.

###  II-C Solving the Polar Problem

Specifically, for our regularization function (11), we show the following
result which allows the polar problem to be solved efficiently, enabling
efficient and guaranteed optimization (14).

######  Theorem 1.

For (14), the polar problem in (16) is

$\displaystyle\Omega_{\theta}^{\circ}(\mathbf{Z})=$
$\displaystyle\max_{\mathbf{d},\mathbf{x},k}\,\mathbf{d}^{\top}\mathbf{Z}\mathbf{x}$
$\displaystyle\mathrm{s.t.}~{}\mathbf{d}^{\top}\mathbf{A}(k)\mathbf{d}\leq
1,\|\mathbf{x}\|^{2}_{F}\leq 1,\ 0\leq k\leq 2.$

where $\mathbf{A}(k)$ is as defined in (13). Further, if we define $k^{*}$ as

$k^{*}=\operatorname*{arg\,max}_{k\in[0,2]}\|\mathbf{A}(k)^{-1/2}\mathbf{Z}\|_{2}\;.$
(18)

Then the optimal values of $\mathbf{d},\mathbf{x},k$ are given as
$\mathbf{d}^{*}=\mathbf{A}(k^{*})^{-1/2}\widebar{\mathbf{d}}$,
$\mathbf{x}^{*}=\widebar{\mathbf{x}}$, and $k^{*}$. Where
$\widebar{\mathbf{d}}$ and $\widebar{\mathbf{x}}$ are the left and right
singular vectors, respectively, associated with the largest singular value of
$\mathbf{A}(k^{*})^{-1/2}\mathbf{Z}$.

The proof of Theorem 1 is detailed in supplementary material VII-B. The result
implies that we can solve the polar by performing a one-dimensional line
search over $k$. Due to the fact that the largest singular value of a matrix
is a Lipschitz continuous function, this line search can be solved efficiently
by a variety of global optimization algorithms. In particular, the above line
search over $k$ is Lipschitz continuous (with respect to $k^{2}$) with a
Lipschitz constant, $L_{k}$, which is bounded by

$\displaystyle L_{k}\\!$
$\displaystyle\leq\\!\begin{Bmatrix}\frac{2}{3\sqrt{3}}\sqrt{\gamma}\|\mathbf{Y}-\mathbf{D}\mathbf{X}^{\top}\|_{2}&\gamma\geq\frac{1}{32}\\\
4\gamma(1+16\gamma)^{-\tfrac{3}{2}}\|\mathbf{Y}-\mathbf{D}\mathbf{X}^{\top}\|_{2}&\gamma<\frac{1}{32}\end{Bmatrix}$
$\displaystyle\qquad\quad\leq\tfrac{2}{3\sqrt{3}}\sqrt{\gamma}\|\mathbf{Y}-\mathbf{D}\mathbf{X}^{\top}\|_{2}\\!$
(19)

This is formally stated and proven proven in supplementary material VII-C as
Theorem 2, where we also give a formal result, in Corollary 1, which
guarantees that the line search over $k$ (and hence the overall polar problem)
can be solved in polynomial time. Building from this key result, next we
discuss an algorithm to provably solve our wave-informed factorization model
in polynomial time.

###  II-D Wave-Informed Factorization Algorithm

Algorithm 1 defines our wave-informed matrix factorization algorithm. The
algorithm has three components that are iterated: (1) perform gradient descent
with a fixed number of modes/columns, $N$, to reach a first order stationary
point for $\mathbf{D}$, $\mathbf{X}$, and $k$, (2) solve the polar problem via
Theorem 1 to obtain $\mathbf{d}^{*}$, $\mathbf{x}^{*}$, and $k^{*}$, and (3)
if the stopping condition is not met, update $\mathbf{D}$, $\mathbf{X}$, and
$k$ by appending the polar solution $(\mathbf{d}^{*},\mathbf{x}^{*},k^{*})$ to
the solution (scaled by an optimal step size $\tau$); otherwise, terminate the
algorithm.

Algorithm 1 Wave-Informed Matrix Factorization

1: Input $\mathbf{D}_{init}$, $\mathbf{X}_{init}$, $\mathbf{k}_{init}$,
stopping tolerance: $\epsilon$

2: Initialize
$\left(\mathbf{D},\mathbf{X},\mathbf{k}\right)\leftarrow\left(\mathbf{D}_{init},\mathbf{X}_{init},\mathbf{k}_{init}\right)$

3: Initialize
$\left(\tilde{\mathbf{D}},\tilde{\mathbf{X}},\tilde{\mathbf{k}}\right)\leftarrow\left(\mathbf{D}_{init},\mathbf{X}_{init},\mathbf{k}_{init}\right)$

4: while $\left|\
\Omega^{\circ}_{\theta}(\tfrac{1}{\lambda}(\mathbf{Y}-\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top}))-1\
\right|\geq\epsilon$ do

5: Perform gradient descent on (14) to obtain
$(\widetilde{\mathbf{D}},\widetilde{\mathbf{X}},\widetilde{\mathbf{k}})$

6: Calculate
$\Omega^{\circ}_{\theta}(\tfrac{1}{\lambda}(\mathbf{Y}-\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top}))$
and obtain $\mathbf{d}^{*},\mathbf{x}^{*},k^{*}$

7: if
$\Omega^{\circ}_{\theta}(\tfrac{1}{\lambda}(\mathbf{Y}-\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top}))-1>\epsilon$
then

8:
$\left(\mathbf{D},\mathbf{X},\mathbf{k}\right)\\!\leftarrow\\!\left(\left[\widetilde{\mathbf{D}},\
\tau\mathbf{d}^{*}\right]\\!,\\!\left[\widetilde{\mathbf{X}},\
\tau\mathbf{x}^{*}\right]\\!,\\!\left[\widetilde{\mathbf{k}}^{\top},\
k^{*}\right]^{\top}\right)\\!$

9: end if

10: end while

####  II-D1 Gradient Descent Update

We begin our algorithm by performing block coordinate descent with respect to
$\mathbf{D}$ and $\mathbf{X}$ along with a block minimization step with
respect to $k$ on the objective in (4), which we note is guaranteed to reach a
first-order stationary point [45]. The number of columns of $\mathbf{D}$ and
$\mathbf{X}$ (denoted by $N$ in (14)) is fixed and does not change during this
step. Specifically, we use the following update equations (for a step size
$\alpha$)

$\displaystyle\mathbf{d}_{i}^{+}$
$\displaystyle=\mathbf{d}_{i}-\alpha\left(\left(\mathbf{D}\mathbf{X}^{\top}-\mathbf{Y}\right)\mathbf{x}_{i}+\lambda\mathbf{d}_{i}\right.$
$\displaystyle\qquad\quad\left.+2\gamma\lambda\left(\mathbf{L}+k_{i}^{2}\mathbf{I}\right)^{2}\mathbf{d}_{i}\right)$
(20) $\displaystyle\mathbf{x}_{i}^{+}$
$\displaystyle=\mathbf{x}_{i}-\alpha\left(\left(\mathbf{D}^{+}\mathbf{X}^{\top}-\mathbf{Y}\right)^{\top}\mathbf{d}_{i}^{+}+\lambda\mathbf{x}_{i}\right)$
(21) $\displaystyle k^{+}_{i}$
$\displaystyle=\sqrt{-\frac{(\mathbf{d}^{+}_{i})^{\top}\mathbf{L}\mathbf{d}^{+}_{i}}{\|\mathbf{d}^{+}_{i}\|_{2}^{2}}}\qquad\textrm{for
each $i=1,\ldots,N$}$ (22)

where $\mathbf{D}^{+}$, $\mathbf{X}^{+}$, and $k_{i}^{+}$ are the updated
values and the $i$ subscript denotes a column of the respective matrix.

####  II-D2 Solve Polar Problem

After finding a first order stationary point
$(\widetilde{\mathbf{D}},\widetilde{\mathbf{X}},\widetilde{\mathbf{k}}$, we
next solve the polar problem according to Theorem 1 with
$\mathbf{Z}=\tfrac{1}{\lambda}(\mathbf{Y}-\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top})$.
In particular, we first solve the one-dimensional line search for the optimal
value of $k$ in (18), where this line search can be solved using many
approaches. For example, the LIPO algorithm given in [46] can provably solve
the line search in polynomial time. Once the optimal $k^{*}$ is obtained, the
values of $\mathbf{d}^{*}$ and $\mathbf{x}^{*}$ are then computed directly
from the eigen-decomposition of $\mathbf{A}(k^{*})$ as described in Theorem 1.

####  II-D3 Stopping Condition and Growing the Factorization

After solving the polar problem as described above, the value of the polar
gives a test for checking global optimality, and if this test is not met then
we can escape the first order stationary point by appending the polar solution
to our factorization. Specifically, once we evaluate the polar
$\Omega^{\circ}_{\theta}(\tfrac{1}{\lambda}(\mathbf{Y}-\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top}))$,
then in [38] (Prop. 4) it is shown that the distance to the global optimum (in
objective value) is directly proportional to the value of the polar value
minus 1 (and that the value of the polar is always $\geq 1$), so choosing to
stop the algorithm when the value of polar takes a value $\leq 1+\epsilon$
also guarantees optimality to within $\mathcal{O}(\epsilon)$. If the value of
the polar does not meet our stopping condition, then we continue the
optimization by appending the polar solution to our factorization after an
optimal scaling. Specifically, let

$\left(\mathbf{D}_{\tau},\mathbf{X}_{\tau},\mathbf{k}_{\tau}\right)\\!=\\!\left(\left[\widetilde{\mathbf{D}},\
\tau\mathbf{d}^{*}\right]\\!,\\!\left[\widetilde{\mathbf{X}},\
\tau\mathbf{x}^{*}\right]\\!,\\!\left[\widetilde{\mathbf{k}}^{\top},\
k^{*}\right]^{\top}\\!\right)$ (23)

be our optimization variables with the new polar solution appended to the
factorization scaled by $\tau$, then we find the optimal value of $\tau$ to
minimize the original objective function:

$\displaystyle\min_{\begin{subarray}{c}\tau\end{subarray}}\\!\tfrac{1}{2}\|\mathbf{Y}-\mathbf{D}_{\tau}\mathbf{X}_{\tau}^{\top}\|_{F}^{2}\\!+\\!\tfrac{\lambda}{2}\\!\|\mathbf{X}_{\tau}\|_{F}^{2}\\!+\\!\tfrac{\lambda}{2}\\!\sum_{i=1}^{N}\mathbf{D}_{\tau}^{\top}\mathbf{A}(k_{i,\tau})\mathbf{D}_{\tau}$
(24)

It is shown in supplementary material VII-D that the optimal $\tau$ is given
by

$\displaystyle\tau=\frac{\sqrt{(\mathbf{d}^{*})^{\top}\left(\mathbf{Y}-\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top}\right)\mathbf{x}^{*}-\lambda}}{\|\mathbf{d}^{*}\|_{2}\|\mathbf{x}^{*}\|_{2}}$
(25)

After updating the variables as in (24) with the optimal choice of $\tau$, we
then return to gradient descent updates until reaching a first order
stationary point.

Taken together, our proposed algorithm results in a polynomial-time algorithm
for solving our main objective (14). We do not give a formal result as the
algorithm given in Algorithm 1 is simply one choice of optimization algorithm
which is possible due to the polynomial time solution for the polar problem
that we give in Theorem 1 and the overall convergence of the algorithm is
largely known from prior work, provided the polar problem can be solved. For
example, the authors of [45] prove that the block coordinate gradient descent
portion of the algorithm will converge to a stationary point at a rate of
$\mathcal{O}(1/t)$ (with $t$ being the number of iterations), which can be
improved to $\mathcal{O}(1/t^{2})$ with Nesterov acceleration. Likewise, the
authors of [39] show that a Frank-Wolfe/conjugate-gradient step (i.e., when we
escape poor stationary points by solving the polar problem) will reduce the
objective with a linear convergence rate. From this it is straight-forward to
combine these arguments to give an overall guarantee of convergence, but we do
not do so here to limit the scope and length of the manuscript.

###  II-E Signal Processing Interpretation of the Regularizer

We now discuss an interesting interpretation of the above algorithm from a
signal processing perspective. Specifically, if we let
$\mathbf{L}=\Gamma\Lambda\Gamma^{\top}$ denote a singular value decomposition
of $\mathbf{L}$, then note that when identifying the optimal $k$ value in the
polar program, we solve for

$\operatorname*{arg\,max}_{k\in[0,2]}\|\Gamma(\mathbf{I}+\gamma(k^{2}\mathbf{I}+\Lambda)^{2})^{-1/2}\Gamma^{\top}\mathbf{Z}\|_{2}\;.$
(26)

This optimization has an intuitive interpretation from signal processing.
Given that $\Gamma$ contains the eigenvectors of a Toeplitz matrix, those
eigenvectors have spectral qualities similar to the discrete Fourier transform
(the eigenvectors of the related circulant matrix would be the discrete
Fourier transform [47, 48]). As a result, $\Gamma^{\top}$ transforms the data
$\mathbf{Z}$ into a spectral-like domain and $\Gamma$ returns the data back to
the original domain. Since the other terms are all diagonal matrices, they
represent element-wise multiplication across the data in the spectral domain.
This is equivalent to a filtering operation, with filter coefficients given by
the diagonal entries of
$(\mathbf{I}+\gamma(\bar{k}\mathbf{I}+\Lambda)^{2})^{-1/2}$. Figure 1 shows
two such filter responses at $k^{2}=2.5,\gamma=1000$ and
$k^{2}=1.5,\gamma=10000$. Observe that, an increased value of $\gamma$ has
reduced the bandwidth of the filter.

Figure 1: Butterworth Filter on the spectral-like domain

Furthermore, recall that the transfer function of a 1st-order Butterworth
filter is given by:

$T(\omega)=\frac{1}{\sqrt{1+\gamma(\omega_{0}+\omega)^{2}}}$ (27)

where $\omega_{0}$ is the center frequency of the passband of the filter and
$1/\sqrt{\gamma}$ corresponds to the filter’s $-3$dB cut-off frequency.
Comparing this to the filter coefficients from (26), we note that the filter
coefficients are identical to those of the 1st-order band-pass Butterworth
filter, where $\Lambda$ corresponds to the angular frequencies.

As a result, we can consider this optimization as determining the optimal
filter center frequency $\omega_{0}=k^{2}$ with fixed bandwidth
$(1/\sqrt{\gamma})$ that retains the maximum amount of signal power from
$\mathbf{Z}$. Likewise, the choice of the $\gamma$ hyperparameter sets the
bandwidth of the filter. As $\gamma\rightarrow\infty$, the filter bandwidth
approaches 0 and thereby restricts us to a single-frequency (i.e., Fourier)
solution. Furthermore, we can provide a recommended lower bound for $\gamma$
according to $\gamma>1/k_{bw}^{2}$, where $k_{bw}$ is the bandwidth of the
signal within this spectral-like domain. We would like to emphasize that this
work introduces a unique theme within the works of scientific machine
learning. Through means of wave-informed learning, we obtain a new method to
mathematically derive the well-known _digital Butterworth filter_.

##  III Experimental Evaluation

(a) (b) (c) (d)

Figure 2: (a) Homogeneously vibrating medium marked with different colors at
three different time instances for a short spatial interval at the start. (b)
Non-Homogeneously (exponentially decaying over space) vibration data marked
with different colors at three different time instances for a short spatial
interval at the start. (c) Traveling wave data marked with different colors at
three different time instances for a short spatial interval at the start. (d)
Data of two mediums joined end-to-end and vibrating together marked with
different colors at different time instances for a short spatio-temporal
interval around the joint.

We demonstrate the results from wave-informed matrix factorization on four
simulated vibration datasets, each characterized by increasing complexity: (1)
a homogeneous vibration with fixed boundary conditions, (2) an non-homogeneous
vibration with fixed boundary conditions and spatial amplitude decay, (3) a
plane wave traveling with temporal decay, and (4) a vibration in a segmented
medium, where each segment has an independent wavenumbers and amplitude. In
all four datasets, one dimension of the matrix represents space and the other
dimension represents time, as is common in modal analysis [49, 50]. Multiple
modes share the same spatial locations. Fig. 2 illustrates these four
datasets.

For each method and dataset, we study five different signal-to-noise ratio
(SNR) levels. Each signal is defined by

$\displaystyle\widehat{y}(\ell,t)=y(\ell,t)+\eta(\ell,t)\;,$ (28)

where $y(\ell,t)$ represents our data and $\eta(\ell,t)$ represents our noise.
The noise is assumed to be additive white Gaussian noise across both space and
time. The SNR is defined by

$\displaystyle\textrm{SNR}=\frac{\sum_{\ell,t}|y(\ell,t)|^{2}}{\sum_{\ell,t}|\eta(\ell,t)|^{2}}\;.$
(29)

The SNR is assumed to be unknown to each of the algorithms.

The following subsections describe how each of our datasets are computed.
Table I lists the parameters for each of these scenarios. Note that each case
can have potentially an infinite number of modes $N$. In practice, these modes
usually attenuate as the associated frequency / wavenumber increases. To
simplify this, we simulate a finite number of modes for each dataset. For the
segmented data, $k_{n}$ represents the first $8$ theoretical modes (four in
each segment) for vibration across two connected segments.

To solve the polar problem and estimate $k^{*}$ in Algorithm 1, many
approaches can be taken. In this paper, we solve for $k^{*}$ by first coarsely
searching for a maximum among coarsely uniformly placed points. For a
predetermined region around the maximum, we perform a fine search by taking
the maximum over finely uniformly placed points in that region

TABLE I:  Table of parameters used to generate data, with additive noise $\eta_{n}\sim\mathcal{N}(0,1)$ and $u_{n}\sim\mathcal{U}[0,1]$. The random variables are independently sampled for each parameter and $1\leq n\leq N$. Data | $(\Delta\ell,L)$ | $(\Delta t,T)$ | $\alpha_{n}$ | $\beta_{n}$ | $k_{n}$ | $\omega_{n}$ | N | $\delta$
---|---|---|---|---|---|---|---|---
Homogeneous Vibration | $(0.0901,1)$ | $(0.0005,2)$ | $\frac{n}{2}$ | $0$ | $n\pi$ | $106k_{n}$ | 6 | 100
Non-homogeneous Vibration | $(0.901,10)$ | $(0.0005,2)$ | 0 | $\frac{4+n}{10}+u_{n}$ | $n\pi$ | $106k_{n}$ | 6 | 1
Travelling Plane Wave | $(0.901,10)$ | $(0.0005,2)$ | $\frac{n}{2}$ | 0 | $n\pi$ | $106k_{n}$ | 6 | 1
Segmented Vibration | $(0.01,2)$ | $(0.02,10)$ | - | - | $2\pi\pm\arctan\left(\sqrt{\frac{1}{3}\left(13\pm 4\sqrt{10}\right)}\right)$ | $3k_{n}$ | $8$ | 1

###  III-A Homogeneous Vibration Data with Fixed Boundaries

Our first dataset emulates vibrations in a medium with fixed (Dirichlet)
boundary conditions such that $y(0,t)=0$ and $y(L,t)=0$, where $L$ is the
length of the medium. The sinusoids reflect from each end and create a
standing wave in space. Such data is defined by

$\displaystyle
y(\ell,t)\\!=\\!\sum_{n=1}^{N}e^{-\alpha_{n}t}\sin(nk\ell)\left(a_{n}\sin(n\omega
t)\\!+\\!b_{n}\cos(n\omega t)\right)\,.$ (30)

This type of data is commonly found in vibrations and modal analysis problems
[51, 52]. For example, the modes and their parameters can be used to determine
material or structural characteristics [8, 9, 10]. Apart from the added noise,
we assume a damping over time to introduce a non-ideality. For this
illustration we choose $L=1$.

###  III-B Inhomogeneous Vibration Data With Spatial Decay

The previous dataset assumed spatial homogeneity. That is, the behavior of the
medium (i.e., the wavenumber and amplitude) did not vary across space $\ell$.
This assumption may not always be true. As a result, we would like to extract
modes with inhomogenieties. We consider the same fixed boundary data but with
a spatial decay such that

$\displaystyle
y(\ell,t)=\sum_{n=1}^{N}e^{-\beta_{n}\ell}\sin(k_{n}\ell)\left(a_{n}\sin(n\omega
t)+b_{n}\cos(n\omega t)\right)\,.$ (31)

In this expression, we have introduced a new parameter $\beta_{n}$ that causes
the amplitude to decay over space, similar to an evanescent wave in electro-
magnetics and is usually meant to persist over longer regions of space. This
example is meant to explore if each algorithm addresses amplitude non-
homogeneities. In contrast to the previous case, where we considered measuring
vibrations along a string of length 1, we consider a scenario of measuring
waves for a longer length. In this illustration, we fix a length of 10 for
measuring waves decaying exponentially over space. As an illustration for the
advantage of length, the spectrum (computed by multiplying data with the
transpose of the eigenvectors $\Gamma$) of three exponentially decaying
sinusoids (for length 1 and length 10) corresponding to $k_{1}=\pi$,
$k_{2}=2\pi$, and $k_{3}=3\pi$ are shown in Figure 3(a) and Figure 3(b),
respectively. In each plot, the peaks correspond to the same three values of
$k$, but due to the 10 times increases in length, Figure 3(b) exhibits a 10
times improvement in resolution in the wavenumber domain. In the presence of
spectral spread (e.g., caused by a spatial decay), this can improve our
ability to distinguish each mode. The parameters generating the data are given
in Table I.

(a) (b)

Figure 3:  Spectra (with respect to Eigen-basis of $\mathbf{L}$) of decaying
exponentials $e^{-\beta_{1}\ell}\sin(\pi\ell)$,
$e^{-\beta_{2}\ell}\sin(2\pi\ell)$ and $e^{-\beta_{3}}\sin(3\pi\ell)$
(indicated using blue, red and yellow colors respectively) for (a)
$\ell\in[0,1]$, since here L=1 (b) $\ell\in[0,10]$, since here L=10

###  III-C Travelling Plane Wave Data

We previously considered vibration-like modes with fixed boundaries. This
dataset considers sinusoids of uniformly spaced wavenumbers travelling in time
and space without explicit boundary conditions, defined by

$y(\ell,t)=\sum_{n=1}^{N}\,e^{-\alpha_{n}t}\sin(k_{n}\ell+\omega_{n}t)\;.$
(32)

This expression resembles a travelling plane wave rather than a vibrational
wave. These types of data are commonly found in beamforming problems for RADAR
[11, 53], ultrasound [54, 13], audio processing [55] and other fields using
spatial signal processing. In beamforming, the parameter $k_{n}$ is usually
proportional to the direction of arrival. In ultrasound and geophysics,
$k_{n}$ may represent a transverse or longitudinal modes of vibration [56]. We
also have the same assumption of measuring the traveling wave data for a
longer length (here 10 times the homogeneous case). Note that time and space
are not immediately separable in this equation, as observed in the other
datasets. However, (32) can be expressed as

$\displaystyle y(\ell,t)$
$\displaystyle=\sum_{n=1}^{N}\,e^{-\alpha_{n}t}\sin(k_{n}\ell-\phi_{n}+\phi_{n}+\omega_{n}t)$
$\displaystyle=\sum_{n=1}^{N}\,e^{-\alpha_{n}t}\left[\sin(k_{n}\ell-\phi_{n})\cos(\omega_{n}t+\phi_{n})\right.$
$\displaystyle\qquad\qquad\qquad+\left.\cos(k_{n}\ell-\phi_{n})\sin(\omega_{n}t+\phi_{n})\right]$
(33)

for an any arbitrarily chosen $\phi_{n}$. Hence, a travelling plane wave mode
can be expressed as the sum of two standing wave modes, similar to our second
dataset, with 90 degree phase differences. In our results, we assume
$\phi_{n}=0$ so that errors in Table III will be small when the extracted
components match $\phi_{n}=0$ while errors in Table III will be invariant to
$\phi_{n}$ since the magnitude of the Fourier transform ignores phase delay.

###  III-D Multi-Segment Data

Our final dataset considers a problem with both amplitude and wavenumber
spatial inhomogeneities. This occurs when there are media adjacent to each
other. As a signal travels from one medium to the next, the amplitude varies
according the the interface’s reflection and transmission coefficients and the
wavenumber varies according to Snell’s law. We generate vibration data in two
media connected end-to-end that is continuous and once differentiable at the
point where the two different they connect (see [56] for related theory). We
assume the two mediums to have a density ratio of 9 to 1, respectively. The
velocity of the wave traveling in the string is inversely proportional to the
square root of density, thus the velocities must be in the ratio 3 to 1. The
mathematical expression generating the data is:

$\displaystyle
y(\ell,t)=\begin{cases}\sum_{i=1}^{N}\sin(3k_{i}\ell)\sin(\omega_{i}(t/c)),&0\leq\ell<1\\\
\sum_{i=1}^{N}\sin(k_{i}\ell)(\ell)\sin(\omega_{i}(t/3c)),&1<\ell\leq
2\;.\end{cases}$ (34)

where $c$ is the velocity of propagation in the first medium. The wavenumbers
for this problem take a specific solution based on the condition, shown in the
Table I, assuming fixed, zero-valued boundaries at $\ell=0$ and $\ell=2$ and
the continuity condition at $\ell=1$.

###  III-E Choice of Regularization Parameters

To provide a fair and informative comparison between each method,we assume
oracle knowledge of $N$ throughout our results. As previously mentioned,
deviation from $1$ in the polar problem in wave informed wave factorization
can be used as a metric for monitoring convergence and as a stopping condition
for the algorithm (and by extension the final choice of $N$), but the
appropriate threshold is unclear in the presence of noise and a similar
convergence property does not exist for all of our comparison methods. While
there are a number of statistical methods for estimating $N$ based on singular
value decomposition that could be applied to every method [57, 58], we want
our performance to reflect choice of algorithms rather than errors in the
estimation of $N$. Hence using oracle knowledge of $N$ to stop the algorithm
as appropriate.

For wave-informed matrix factorization, we select $\lambda$ to be
$3/4\sigma(N)$, where $\sigma(N)$ is the N-th singular value of our data
matrix $\mathbf{Y}$. This choice of $\lambda$ is intuitively chosen based on
the interpretation of low rank factorization as a soft thresholding operation
over the singular values [37, 38]. The $3/4$ factor was obtained empirically
based on our examples. However, reducing the factor is reasonable given the
additional constraint provided by the Helmholtz equation.

We select $\gamma$ in wave-informed matrix factorization based on our
intuition of $1/\sqrt{\gamma}$ being bandwidth of a filter. In general, the
bandwidth of exponentially decaying sinusoids discussed here is small, so
$\gamma$ is expected to be relatively large. Since our modes are constrained
to a finite length in space, we know the bandwidth of a pure sinusoid will be
one sample of the discrete Fourier transform. Therefore, we choose our
bandwidth based on this assumption so that

$\displaystyle\gamma=\delta\left(\frac{M}{\pi}\right)^{2},\;$ (35)

where $M$ is the number of samples in the spatial domain. While this bandwidth
is slightly smaller than some of the data due to spatial variations, this
difference should be corrected for in the gradient descent step of wave-
informed matrix factorization. An analytical way to motivate the same would be
to observe the filter coefficients.

$\displaystyle\cfrac{1}{\sqrt{1+\gamma(k^{2}+\boldsymbol{\Lambda}_{ii})^{2}}}$

is the $i$-th filter coefficient. Assuming pure Dirichlet boundary conditions,
we have from [43] that,

$\displaystyle\boldsymbol{\Lambda}_{ii}=-4\sin^{2}\left(\cfrac{\pi
i}{2(M+1)}\right)$

For $M$ large enough compared to $i$, we have that,

$\displaystyle\boldsymbol{\Lambda}_{ii}=-4\sin^{2}\left(\cfrac{\pi
i}{2(M+1)}\right)\approx-\cfrac{\pi^{2}i^{2}}{(M+1)^{2}}$

Again, for large enough $M$, we can approximate the filter coefficients as
below by substituting $\gamma$ from (35),

$\displaystyle\cfrac{1}{\sqrt{1+\delta\cfrac{M^{2}}{\pi^{2}}\left(k^{2}-\cfrac{\pi^{2}i^{2}}{(M+1)^{2}}\right)^{2}}}\approx\cfrac{1}{\sqrt{1+\delta\left(\cfrac{M^{2}k^{2}}{\pi^{2}}-i^{2}\right)^{2}}}$

where $i$ represents the actual (not angular) wavenumber scale. For instance,
when $\delta=1$, the bandpass Butterworth filter has a -3 dB cut-off
wavenumber of 1 m-1 in the actual wavenumber scale. For larger values of
$\delta$, the -3dB cut-off wavenumber reduces to $1/\delta$ m-1 and is useful
to pick isolated wavenumbers. The value of $\delta$ used for each case is
mentioned in Table I. Observe that homogeneous vibrations requires a larger
value of $\delta(=100)$ to pick isolated wavenumbers whereas others still work
with a value of $\delta=1$.

##  IV Results & Discussions

For each dataset, we compare our approach with five other non-parametric
matrix factorization / blind source separation methods. The algorithms we
compare with include independent component analysis (ICA) [59], dynamic mode
decomposition (DMD) [60], multidimensional empirical mode decomposition (EMD)
[61, 62, 63], principal component analysis (PCA) [64], and wave-informed
K-singular value decomposition (WIKSVD) [31]. We choose these method because
they are widely used in modal analysis throughout the literature or represent
a precursor to wave-informed matrix factorization (WIMF). Each method places
different assumptions on the extracted modes. Briefly, ICA extracts components
components that are statistically independent, PCA extracts components that
are orthogonal to each other, DMD extracts components with a set frequency and
decay based on estimating a state transition matrix, EMD extracts components
by iteratively interpolating between maxima and minima, and WIKSVD extracts
components that assume spatial components satisfy the wave equation and are
sparse in the frequency domain. By comparison, WIMF extracts spatial
components that satisfy the wave equation and minimizes the number of
extracted modes. Based on the filter interpretation of WIMF, it has
characteristics similar to DMD since the filter’s center frequency and
bandwidth identifies the frequency and possibly decay of the extracted modes.

###  IV-A General Performance

In Table III, we show the mean squared errors for each of the recovered modes
obtained by each method,

$\displaystyle\textrm{MSE}=\frac{1}{N}\sum_{n=1}^{N}\left\|\frac{\mathbf{d}_{n}}{\|\mathbf{d}_{n}\|_{2}}-\frac{c_{n}\mathbf{d}_{n}^{(\textrm{true})}}{\|\mathbf{d}_{n}^{(\textrm{true})}\|_{2}}\right\|^{2}_{2}$
(36)

where $\mathbf{d}_{n}$ is extracted mode and
$\mathbf{d}_{n}^{(\textrm{true})}$ is the true mode. In our analysis,
$\mathbf{d}_{n}^{(\textrm{true})}$ is chosen (without replacement) to be the
mode with the maximal correlation (in absolute value) with $\mathbf{d}_{n}$
and
$c_{n}=\textrm{sign}(\mathbf{d}_{n}^{\top}\mathbf{d}^{\textrm{(true)}}_{n})$.
When computing the mean squared error, we normalize each mode since the
amplitude is not uniquely determined in matrix decompositions.

Table III shows the mean squared error between the Fourier magnitudes true and
extracted modes, using the same definition of (36) but replacing
$\mathbf{d}_{n}$ and $\mathbf{d}_{n}^{(\textrm{true})}$ with the magnitude of
their Fourier transforms, respectively. These results provide a different
perspective, identifying if the methods obtain the correct frequencies
ignoring the effect of phase shifts.

Fig. 5 further illustrates examples of two modes extracted by each method for
each dataset at an SNR of $\infty$. Across all of the datasets, ICA has the
overall poorest performance, followed by EMD, PCA, DMD, WIKSVD, and WIMF. We
hypothesize ICA performs poorly since wave modes are not random signals, and
while the modes are orthogonal, they are not statistically independent. EMD
performs poorly since mixtures of sinusoids where the amplitudes and/or
frequencies are too similar cannot be recovered by the algorithm [65]. DMD and
WIKSVD performs most poorly in high noise scenarios and with spatially
segmented wavenumbers. DMD performs poorly at high noise because the
estimation of the state transition matrices does not consider noise. Though a
similar wave-informed regularizer is introduced for WIKSVD, it does not
extract the basis as well as WIMF we believe due to the lack of filtering
behavior in WIKSVD. This filtering effect particularly improves WIMF over
WIKSVD for high noise regimes. The following subsections study the performance
with each dataset in greater detail.

TABLE II:  Mean squared errors (times $10^{3}$) between actual basis
generating the data for each case (averaged over 20 Monte-Carlo experiments)
and the matrix $\mathbf{D}$ obtained from the algorithms (WIMF), (WIKSVD),
(ICA), (DMD), (EMD) and (PCA). Highlight values show the best performance.

Data | SNR [dB] | WIKSVD | ICA | DMD | EMD | PCA | WIMF
---|---|---|---|---|---|---|---
Homogeneous | $-11.84$ | p m $0.38$ | $0.77$ | p m $0.13$ | $1826.96$ | p m $14.19$ | $233.9$ | p m $18.05$ | $1173.39$ | p m $159.52$ | $61.78$ | p m $8.21$ | $0.69$ | p m $0.45$
Homogeneous | $-5.85$ | p m $0.55$ | $0.18$ | p m $0.06$ | $1823.34$ | p m $11.24$ | $63.94$ | p m $7.98$ | $1132.57$ | p m $130.47$ | $13.17$ | p m $2.1$ | $0.21$ | p m $0.08$
Homogeneous | $-1.93$ | p m $0.31$ | $0.08$ | p m $0.01$ | $1825.9$ | p m $11.85$ | $26.89$ | p m $2.2$ | $1206.98$ | p m $151.27$ | $5.25$ | p m $0.65$ | $0.14$ | p m $0.1$
Homogeneous | $2.28$ | p m $0.36$ | $0.04$ | p m $0.01$ | $1824.56$ | p m $10.09$ | $10.5$ | p m $1.06$ | $1243$ | p m $133.39$ | $2.03$ | p m $0.27$ | $0.09$ | p m $0.05$
Homogeneous | $\infty$ |  | $0.018$ | p m $0.001$ | $1155.081$ | p m $96.941$ | $0.172$ | p m $0.022$ | $1443.895$ | p m $160.807$ | $0.063$ | p m $0.024$ | $0.057$ | p m $0.009$
Non-homogeneous | $-13.09$ | p m $0.59$ | $7.56$ | p m $1.1$ | $1917.6$ | p m $12.35$ | $252.93$ | p m $26.97$ | $1402.98$ | p m $181.83$ | $94.83$ | p m $60.72$ | $3.91$ | p m $0.57$
Non-homogeneous | $-9.41$ | p m $0.73$ | $2.99$ | p m $0.38$ | $1920.05$ | p m $6.36$ | $133.43$ | p m $17.47$ | $1414.42$ | p m $189.53$ | $48.07$ | p m $48.87$ | $3.39$ | p m $0.45$
Non-homogeneous | $-3.38$ | p m $0.68$ | $0.76$ | p m $0.1$ | $1922.28$ | p m $7.17$ | $38.64$ | p m $5.61$ | $1390.92$ | p m $138.97$ | $12.54$ | p m $19$ | $3.12$ | p m $0.36$
Non-homogeneous | $2.69$ | p m $0.56$ | $0.21$ | p m $0.03$ | $1922.94$ | p m $5.0$ | $10.42$ | p m $1.75$ | $1381.22$ | p m $137.57$ | $12.59$ | p m $16.24$ | $3.14$ | p m $0.45$
Non-homogeneous | $\infty$ |  | $0.039$ | p m $0.005$ | $1956.643$ | p m $67.1$ | $0.000\,041$ | p m $0.000\,006$ | $1257.568$ | p m $162.273$ | $25.135$ | p m $33.957$ | $3.042$ | p m $0.37$
Travelling Wave | $-11.73$ | p m $0.627$ | $320.05$ | p m $106.177$ | $1913.03$ | p m $13.517$ | $1170.23$ | p m $57.756$ | $1415.61$ | p m $129.838$ | $334.24$ | p m $82.329$ | $1129.04$ | p m $235.417$
Travelling Wave | $-5.9$ | p m $0.543$ | $225.47$ | p m $78.929$ | $1927.56$ | p m $9.552$ | $1108.11$ | p m $32.655$ | $1346.61$ | p m $115.312$ | $212.16$ | p m $89.363$ | $40.24$ | p m $1.464$
Travelling Wave | $-1.84$ | p m $0.565$ | $199.87$ | p m $49.975$ | $1930.35$ | p m $13.662$ | $1084.06$ | p m $38.305$ | $1267.77$ | p m $110.204$ | $241.56$ | p m $67.66$ | $38.29$ | p m $0.563$
Travelling Wave | $1.83$ | p m $0.834$ | $235.93$ | p m $52.036$ | $1932.83$ | p m $18.775$ | $1094.88$ | p m $29.441$ | $1390.48$ | p m $163.535$ | $280.82$ | p m $62.795$ | $36.66$ | p m $0.408$
Travelling Wave | $\infty$ |  | $409.45$ | p m $13.473$ | $1189.03$ | p m $43.279$ | $1115.94$ | p m $1.251$ | $1548.12$ | p m $75.384$ | $342.49$ | p m $5.63$ | $35.82$ | p m $0.271$
Segmented | $-11.56$ | p m $0.0004$ | $1588.98$ | p m $38.1574$ | $1792.13$ | p m $25.247$ | $1515.05$ | p m $84.4314$ | $1264.24$ | p m $116.8219$ | $1526.52$ | p m $77.3674$ | $891.77$ | p m $117.3127$
Segmented | $-9.07$ | p m $0.0244$ | $1594.08$ | p m $35.9512$ | $1794.89$ | p m $29.2853$ | $1532.5$ | p m $78.7935$ | $1260.97$ | p m $131.089$ | $1491.62$ | p m $60.1776$ | $913.24$ | p m $117.485$
Segmented | $-3.05$ | p m $0.0198$ | $1638.14$ | p m $16.267$ | $1792.46$ | p m $27.3304$ | $1520.51$ | p m $68.7968$ | $1191.45$ | p m $127.6302$ | $1480.75$ | p m $34.9955$ | $767.29$ | p m $20.297$
Segmented | $5.48$ | p m $0.0187$ | $1664.01$ | p m $46.2336$ | $1786.2$ | p m $28.6862$ | $1483.33$ | p m $24.9771$ | $1158.02$ | p m $133.9853$ | $1470.5$ | p m $40.8251$ | $616.28$ | p m $4.7015$
Segmented | $\infty$ |  | $889.88$ | p m $0.001$ | $1840$ | p m $12.7936$ | $1006.11$ | p m $0.001$ | $1019.05$ | p m $0.001$ | $1516.11$ | p m $0.001$ | $852.02$ | p m $0.1228$

TABLE III:  Mean squared errors (times $10^{3}$) between absolute value of
Fourier transform of actual basis generating the data for each case (averaged
over 20 Monte-Carlo experiments) and the matrix $\mathbf{D}$ obtained from the
algorithms (WIMF), (WIKSVD), (ICA), (DMD), (EMD) and (PCA). Highlighted values
show the best performance.

Data | SNR [dB] | WIKSVD | ICA | DMD | EMD | PCA | WIMF
---|---|---|---|---|---|---|---
Homogeneous | $-11.84$ | p m $0.38$ | $0.51$ | p m $0.12$ | $1731.99$ | p m $23.89$ | $229.51$ | p m $18.19$ | $891.27$ | p m $151.17$ | $59.92$ | p m $7.87$ | $0.23$ | p m $0.14$
Homogeneous | $-5.85$ | p m $0.55$ | $0.11$ | p m $0.05$ | $1719.27$ | p m $21.25$ | $62.51$ | p m $7.84$ | $838.18$ | p m $150.21$ | $12.6$ | p m $2.05$ | $0.06$ | p m $0.03$
Homogeneous | $-1.93$ | p m $0.31$ | $0.04$ | p m $0.01$ | $1725.65$ | p m $24.87$ | $26.12$ | p m $2.17$ | $733.01$ | p m $104.9$ | $4.88$ | p m $0.59$ | $0.04$ | p m $0.02$
Homogeneous | $2.28$ | p m $0.36$ | $0.02$ | p m $0.01$ | $1723.94$ | p m $25.25$ | $10.13$ | p m $1.02$ | $768.8$ | p m $86.18$ | $1.83$ | p m $0.21$ | $0.03$ | p m $0.01$
Homogeneous | $\infty$ |  | $0.011$ | p m $0.001$ | $782.744$ | p m $110.157$ | $0.152$ | p m $0.019$ | $1020.455$ | p m $111.294$ | $0.016$ | p m $0.006$ | $0.017$ | p m $0.002$
Non-homogeneous | $-13.09$ | p m $0.59$ | $4.25$ | p m $0.63$ | $1767.38$ | p m $20.51$ | $238.68$ | p m $24.97$ | $958.78$ | p m $122.03$ | $78.86$ | p m $52.51$ | $1.8$ | p m $0.42$
Non-homogeneous | $-9.41$ | p m $0.73$ | $1.5$ | p m $0.19$ | $1771.97$ | p m $15.62$ | $125.96$ | p m $16.54$ | $944.35$ | p m $135.65$ | $38.37$ | p m $39.95$ | $1.47$ | p m $0.31$
Non-homogeneous | $-3.38$ | p m $0.68$ | $0.32$ | p m $0.04$ | $1779.36$ | p m $15.44$ | $35.97$ | p m $5.36$ | $918.64$ | p m $158.62$ | $8.67$ | p m $14.17$ | $1.32$ | p m $0.26$
Non-homogeneous | $2.69$ | p m $0.56$ | $0.08$ | p m $0.01$ | $1771.68$ | p m $15.31$ | $9.45$ | p m $1.68$ | $869.36$ | p m $151.3$ | $8.36$ | p m $12.72$ | $1.34$ | p m $0.3$
Non-homogeneous | $\infty$ |  | $0.017$ | p m $0.002$ | $1753.32$ | p m $79.811$ | $0.000\,006$ | p m $0.000\,001$ | $948.256$ | p m $116.393$ | $1.274$ | p m $0.254$ | $18.944$ | p m $28.261$
Travelling Wave | $-11.73$ | p m $0.627$ | $102.57$ | p m $62.53$ | $1878.85$ | p m $14.601$ | $668.07$ | p m $23.091$ | $1113.04$ | p m $74.405$ | $147.76$ | p m $38.868$ | $628.26$ | p m $187.239$
Travelling Wave | $-5.9$ | p m $0.543$ | $8.91$ | p m $1.364$ | $1893.01$ | p m $8.316$ | $597.02$ | p m $2.652$ | $1040.68$ | p m $64.72$ | $26.04$ | p m $4.451$ | $33.99$ | p m $0.81$
Travelling Wave | $-1.84$ | p m $0.565$ | $5.1$ | p m $0.825$ | $1899.45$ | p m $12.028$ | $588.02$ | p m $1.554$ | $970.5$ | p m $63.797$ | $10.13$ | p m $2.314$ | $32.92$ | p m $0.558$
Travelling Wave | $1.83$ | p m $0.834$ | $3.95$ | p m $1.065$ | $1902.72$ | p m $16.441$ | $585.46$ | p m $0.532$ | $1088.83$ | p m $150.451$ | $4.08$ | p m $0.976$ | $31.6$ | p m $0.329$
Travelling Wave | $\infty$ |  | $1.47$ | p m $0.2$ | $1081.02$ | p m $40.082$ | $585$ | p m $0$ | $1446.47$ | p m $60.782$ | $0.12$ | p m $0.108$ | $31.32$ | p m $0.295$
Segmented | $-11.56$ | p m $0.0004$ | $1037.52$ | p m $36.026$ | $1373.83$ | p m $24.3314$ | $949.2$ | p m $33.3454$ | $724.26$ | p m $80.7343$ | $957.85$ | p m $35.2736$ | $536.7$ | p m $71.3607$
Segmented | $-9.07$ | p m $0.0244$ | $1067.23$ | p m $47.922$ | $1383.48$ | p m $25.9373$ | $907.05$ | p m $27.7681$ | $626.9$ | p m $88.1976$ | $938.54$ | p m $34.8523$ | $502.69$ | p m $98.836$
Segmented | $-3.05$ | p m $0.0198$ | $1249.49$ | p m $79.8907$ | $1382.94$ | p m $25.0862$ | $863.72$ | p m $45.2521$ | $496.23$ | p m $39.492$ | $923.74$ | p m $26.4544$ | $397.71$ | p m $128.3619$
Segmented | $5.48$ | p m $0.0187$ | $1175.7$ | p m $18.6612$ | $1372.06$ | p m $24.5076$ | $973.45$ | p m $37.9708$ | $448.47$ | p m $63.7549$ | $908.77$ | p m $21.7454$ | $428.47$ | p m $6.6073$
Segmented | $\infty$ |  | $623.45$ | p m $0.001$ | $1337.82$ | p m $25.0305$ | $348.7$ | p m $0.001$ | $426.55$ | p m $0.001$ | $887.16$ | p m $0.001$ | $357.87$ | p m $0.1348$

###  IV-B Homogeneous Vibration Data With Fixed Boundaries

Table III and III show that this relatively simple dataset achieves the
overall best performance among the datasets. Wave-informed matrix
factorization, wave-informed K-SVD, and dynamic mode decomposition perform
best, albeit under different conditions. Wave-informed K-SVD and dynamic mode
decomposition perform the best in no-noise conditions. Yet, dynamic mode
decomposition is more sensitive to noise. This is evident in Tables III,III as
with increasing noise the errors indicated in the tables increase. Wave-
informed K-SVD performs best in low to mild noise ($\leq-3$ dB SNR) whilst
wave-informed matrix factorization performs best in heavy noise ($>-3$ dB
SNR).

###  IV-C Inhomogeneous Vibration Data With Spatial Decay

These results follow a similar trend as the homogeneous vibration data. In the
inhomogeneous data, the spatial decay spreads the signal across the wavenumber
domain, resulting in generally less compressible data that is more sensitive
to noise. As a result, we see a reduction in performance for WIMF, WIKSVD, DMD
in Table III and III. In case of infinite SNR, DMD performs best as it is
designed to extract damped sinusoids [66]. However, it remains sensitive to
noise.

PCA and EMD appear unstable with this dataset. This occurs because the
parameter $\beta_{n}$ varies across each method and SNR. This randomness in
conjunction with the incorrect extraction of multiple wavenumbers in a single
component, as previously described for the homogeneous data, is the source for
this error. As a result, EMD and PCA are to be sensitive to the variation in
spatial decay.

Figure 4: (a) Homogeneous Vibration

(b) In-homogeneous Vibration

(c) Traveling Plane Waves

(d) Multiple Segments Figure 5: Two recovered modes (rows) obtained from
(WIMF), (WIKSVD), (PCA), (ICA), (DMD), (EMD) for SNR = $\infty$.

Figure 6: (a) Two columns from the basis obtained from PCA and the absolute
value in the transformed domain (transformed to the eigenbasis of
$\mathbf{L}$) (b) Two columns from the basis obtained from wave-informed
matrix factorization and the absolute value in the transformed domain

###  IV-D Traveling Plane Wave Data

Due to the expression derived in (33), we extract $2N$ modes rather than $N$
modes. This is seen in Fig. 5(c), where the two most dominant modes have the
same frequencies. It is also important to note that the decomposition is not
unique due to the arbitrary phase term $\phi_{n}$. This discrepancy seen in
the difference in PCA’s relatively poor performance in Table III (where
$\phi_{n}$ is assumed to be $0$) and the excellent performance in Table III,
which ignores the phase in its computation. We further observe in Fig. 5(c)
that PCA does not have precise phases of $0$ and $90$ degrees. Similarly, DMD
incorrectly extracts modes with equal phases. This is because DMD extracts a
complex-valued factorization. This reconstructs the data correctly but creates
a large error when compared with the true real-valued modes. Note that the
Fig. 5(c) shows the real part of the DMD result but the result in Table III
and III compute the error with the complex values. In no other datasets does
DMD produce a strong imaginary component.

In addition, WIKSVD appears to have difficulty extracting a mode with a $90$
degree phase based on the errors in Table III and that the second most
prominent mode in Fig. 5(c) has the wrong frequency. This can be explained due
to our choice of Laplacian matrix $\mathbf{L}$, which inherently assumes
fixed, zero-valued boundary conditions. Similarly, we observe that WIMF in
Fig. 5(c) forces zero-valued boundaries, but the second mode immediately
adjusts to have an effective 90 degree phase. This shows that WIMF is more
adaptable than WIKSVD. As a result, WIMF overall achieves the best results,
but has a larger error relative to the previous two datasets.

###  IV-E Multi-Segment Data

Tables III and III show that this dataset is the most difficult of the four.
Under low noise conditions, WIMF achieves the best performance. From Fig.
5(d), we observe the relatively large error is due to wavenumbers information
”leaking” from one segment to the other. In particular, WIKSVD and DMD include
the the modal bevhavior both regions in their components. Hence, these
solution incorrectly include two different wavenumbers in their components.

PCA and EMD achieve decent localization of the modes but fail to extract the
correct frequency content, as evident in Table III. WIMF, EMD, and PCA all
demonstrate some degree of separation and localization, but are not tightly
bound to their spatial region. In particular, WIMF does not achieve stronger
localization because the wave equation constraints are placed on each mode
rather than the complete solution, and the discontinuity between segments does
not satisfy the wave equation for each individual mode.

###  IV-F General discussion on wave-informed matrix factorization

In the performance analysis of this method, we observe that _wave-informed
matrix factorization_ successfully recovers pure sinusoids with zero phase in
the homogeneous vibrations case in Section IV-B and recovers decaying
sinusoids in the inhomogenous vibrations case in Section IV-C. In Section
IV-D, it recovers approximate cosines that start at zero and eventually
converge to the sinusoid shape, when the true result has a phase shift. For
the multi-segment case in Section IV-E, it recovers sinusoidal vibrations that
locally to each region.

We would like to emphasize here that hard constraining the wave equation, i.e.
solving an optimization similar to

$\displaystyle\underset{\mathbf{D},\mathbf{X},N,\mathbf{k}}{\operatorname*{arg\,min}}\frac{1}{2}\|\mathbf{Y}-\mathbf{D}\mathbf{X}^{\top}\|_{F}^{2}+\frac{\lambda}{2}\left(\|\mathbf{D}\|_{F}^{2}+\|\mathbf{X}\|_{F}^{2}\right)$
$\displaystyle\ \ \textnormal{s.t.}\ \
\mathbf{L}\mathbf{D}_{i}=-k_{i}^{2}\mathbf{D}_{i},\;\forall i\in[N]$

would be equivalent to solving, for $\mathbf{D}$, the eigenvectors of
$\mathbf{L}$ that closely span (in the $\ell_{2}$ sense) the column space of
$\mathbf{Y}$. These are pure sinusoids with zero phase. In contrast to this,
we observe that the columns of $\mathbf{D}$ are dictated by both data and the
soft constraint introduced by the wave equation. In particular, we note that
the columns extracted in most of the cases only approximately satisfy the
wave-constraint and do not satisfy the wave equation. This property
distinguishes our algorithm from traditional factorization methods that are
either data dependent or are entirely parameterized by the assumed physics. We
attribute this to the recovery of components such as exponentially decaying
sinusoids (which do not satisfy the wave equation).

##  V Conclusions

We have introduced wave-informed matrix factorization and developed a
framework and an algorithm with provable, global optimally guarantees for the
same. More generally, this work introduced a methodology to enforce linear
homogeneous partial differential equation to influence a matrix factorization
algorithm, and demonstrated this for the case of the time-independent version
of the Helmholtz wave equation. The output from the algorithm was compared
with that of state-of-the-art algorithms for modal and component analysis. We
demonstrated that the wave-informed approach learns representations that are
more physically relevant and practical for the purpose of modal analysis.

Future work will include generalizing this approach to a variety of linear
PDEs beyond the wave equation (especially to adapt the works [52, 49] to our
framework) as well as wave propagation along more than one dimension and
extension to applications in baseline-free anomaly detection for structural
health monitoring [27, 29].

##  VI Acknowledgements

This work is partially supported by NSF EECS-1839704, NSF CISE-1747783, NIH
NIA 1R01AG067396, ARO MURI W911NF-17-1-0304, and NSF-Simons MoDL 2031985.

## References

  * [1] S. Mehrkanoon, “Deep shared representation learning for weather elements forecasting,” _Knowledge-Based Systems_ , vol. 179, pp. 120–128, 2019.
  * [2] M. Ringsquandl, S. Lamparter, R. Lepratti, and P. Kröger, “Knowledge fusion of manufacturing operations data using representation learning,” in _Proc. of the IFIP International Conference on Advances in Production Management Systems_. Springer, 2017, pp. 302–310.
  * [3] N. Zobeiry, J. Reiner, and R. Vaziri, “Theory-guided machine learning for damage characterization of composites,” _Composite Structures_ , vol. 246, p. 112407, 2020.
  * [4] A. Mohamed, G. Hinton, and G. Penn, “Understanding how deep belief networks perform acoustic modelling,” in _Proc. of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2012, pp. 4273–4276.
  * [5] G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. Van Der Laak, B. Van Ginneken, and C. I. Sánchez, “A survey on deep learning in medical image analysis,” _Medical Image Analysis_ , vol. 42, pp. 60–88, 2017.
  * [6] A. Stanbridge and D. Ewins, “Modal testing using a scanning laser doppler vibrometer,” _Mechanical systems and signal processing_ , vol. 13, no. 2, pp. 255–270, 1999.
  * [7] Y. Yang, C. Dorn, T. Mancini, Z. Talken, S. Nagarajaiah, G. Kenyon, C. Farrar, and D. Mascareñas, “Blind identification of full-field vibration modes of output-only structures from uniformly-sampled, possibly temporally-aliased (sub-nyquist), video measurements,” _Journal of Sound and Vibration_ , vol. 390, pp. 232–256, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0022460X16306897
  * [8] D. Kammer, “Sensor placement for on-orbit modal identification and correlation of large space structures,” _Journal of Guidance, Control, and Dynamics_ , vol. 14, pp. 251 – 259, 06 1991.
  * [9] A. Cunha and E. Caetano, “Experimental modal analysis of civil engineering structures,” _Sound and Vibration_ , vol. 40, 06 2006.
  * [10] K. D. Marshall, “Modal analysis of a violin,” _Journal of the Acoustical Society of America_ , vol. 77, pp. 695–709, 1985.
  * [11] J. Capon, “High-resolution frequency-wavenumber spectrum analysis,” _Proceedings of the IEEE_ , vol. 57, no. 8, pp. 1408–1418, 1969.
  * [12] M. Wathelet, B. Guillier, P. Roux, C. Cornou, and M. Ohrnberger, “Rayleigh wave three-component beamforming: signed ellipticity assessment from high-resolution frequency-wavenumber processing of ambient vibration arrays,” _Geophysical Journal International_ , vol. 215, no. 1, pp. 507–523, 07 2018. [Online]. Available: https://doi.org/10.1093/gji/ggy286
  * [13] L. Yu and Z. Tian, “Guided wave phased array beamforming and imaging in composite plates,” _Ultrasonics_ , vol. 68, pp. 43–53, 2016. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0041624X16000160
  * [14] A. Karpatne, G. Atluri, J. H. Faghmous, M. Steinbach, A. Banerjee, A. Ganguly, S. Shekhar, N. Samatova, and V. Kumar, “Theory-guided data science: A new paradigm for scientific discovery from data,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 29, no. 10, pp. 2318–2331, 2017.
  * [15] N. Benjamin Erichson, M. Muehlebach, and M. W. Mahoney, “Physics-informed autoencoders for lyapunov-stable fluid flow prediction,” _arXiv [physics.comp-ph]_ , p. 14, 2019.
  * [16] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” _J. Comput. Phys._ , vol. 378, pp. 686–707, 2019.
  * [17] M. A. Nabian and H. Meidani, “Physics-driven regularization of deep neural networks for enhanced engineering design and analysis,” _Journal of Computing and Information Science in Engineering_ , vol. 20, no. 1, 2019.
  * [18] A. M. Tartakovsky, C. O. Marrero, P. Perdikaris, G. D. Tartakovsky, and D. Barajas-Solano, “Learning parameters and constitutive relationships with physics informed deep neural networks,” _arXiv [math.AP]_ , p. 22, 2018.
  * [19] Y. Long, X. She, and S. Mukhopadhyay, “Hybridnet: Integrating model-based and data-driven learning to predict evolution of dynamical systems,” in _Proc. of the Conference on Robot Learning_ , 2018, p. 10.
  * [20] K. Shukla, P. C. Di Leoni, J. Blackshire, D. Sparkman, and G. E. Karniadakis, “Physics-informed neural network for ultrasound nondestructive quantification of surface breaking cracks,” _arXiv [cs.LG]_ , p. 19, 2020\.
  * [21] A. Karpatne, W. Watkins, J. Read, and V. Kumar, “Physics-guided neural networks (PGNN): An application in lake temperature modeling,” _arXiv [cs.LG]_ , vol. 1710.11431, p. 11, 2017.
  * [22] D. Lei, X. Chen, and J. Zhao, “Opening the black box of deep learning,” _arXiv [cs.LG]_ , vol. 1805.08355, p. 27, 2018.
  * [23] K. R. Mestav, J. Luengo-Rozas, and L. Tong, “State estimation for unobservable distribution systems via deep neural networks,” in _Proc. of the IEEE Power & Energy Society General Meeting_. Portland: IEEE, 2018, pp. 1–5.
  * [24] T.-j. Luo, Y. Fan, L. Chen, G. Guo, and C. Zhou, “EEG signal reconstruction using a generative adversarial network with Wasserstein distance and temporal-spatial-frequency loss,” _Frontiers in Neuroinformatics_ , vol. 14, p. 15, 2020.
  * [25] L. Mosser, O. Dubrule, and M. J. Blunt, “Stochastic seismic waveform inversion using generative adversarial networks as a geological prior,” _Mathematical Geosciences_ , vol. 52, pp. 53–79, 2020.
  * [26] J. B. Harley and J. M. Moura, “Sparse recovery of the multimodal and dispersive characteristics of Lamb waves,” _The Journal of the Acoustical Society of America_ , vol. 133, no. 5, pp. 2732–2745, 2013.
  * [27] K. S. Alguri, J. Melville, and J. B. Harley, “Baseline-free guided wave damage detection with surrogate data and dictionary learning,” _The Journal of the Acoustical Society of America_ , vol. 143, no. 6, pp. 3807–3818, 2018.
  * [28] K. S. Alguri, C. C. Chia, and J. B. Harley, “Model-driven, wavefield baseline subtraction for damage visualization using dictionary learning,” vol. 2, 2017, pp. 2276––2283.
  * [29] ——, “Sim-to-real: Employing ultrasonic guided wave digital surrogates and transfer learning for damage visualization,” _Ultrasonics_ , vol. 111, p. 106338, 2021.
  * [30] K. S. Alguri and J. B. Harley, “Transfer learning of ultrasonic guided waves using autoencoders: A preliminary study,” _AIP Conf. Proc._ , vol. 2102, no. 1, p. 50013, 2019.
  * [31] H. V. Tetali, K. S. Alguri, and J. B. Harley, “Wave Physics Informed Dictionary Learning In One Dimension,” in _Proc. of the IEEE International Workshop on Machine Learning for Signal Processing (MLSP)_. IEEE, 2019, pp. 1–6.
  * [32] H. V. Tetali, J. B. Harley, and B. D. Haeffele, “Wave-informed matrix factorization with global optimality guarantees,” _arXiv preprint arXiv:2107.09144_ , 2021.
  * [33] S. Bhojanapalli, B. Neyshabur, and N. Srebro, “Global optimality of local search for low rank matrix recovery,” in _Advances in Neural Information Processing Systems_ , 2016, pp. 3873–3881.
  * [34] R. Ge, J. D. Lee, and T. Ma, “Matrix completion has no spurious local minimum,” in _Advances in Neural Information Processing Systems_ , 2016, pp. 2973–2981.
  * [35] D. Park, A. Kyrillidis, C. Caramanis, and S. Sanghavi, “Non-square matrix sensing without spurious local minima via the burer-monteiro approach,” _arXiv preprint arXiv:1609.03240_ , p. 14, 2016.
  * [36] R. Ge, C. Jin, and Y. Zheng, “No spurious local minima in nonconvex low rank problems: A unified geometric analysis,” in _Proc. of the International Conference on Machine Learning_ , vol. 70. JMLR. org, 2017, pp. 1233–1242.
  * [37] B. Haeffele, E. Young, and R. Vidal, “Structured low-rank matrix factorization: Optimality, algorithm, and applications to image processing,” in _Proc. of the International Conference on Machine Learning_ , 2014, pp. 2007–2015.
  * [38] B. D. Haeffele and R. Vidal, “Structured low-rank matrix factorization: Global optimality, algorithms, and applications,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 42, no. 6, 2019.
  * [39] F. Bach, “Convex relaxations of structured matrix factorizations,” _arXiv preprint arXiv:1309.3117_ , p. 33, 2013.
  * [40] E. Figueiredo and J. Brownjohn, “Three decades of statistical pattern recognition paradigm for shm of bridges,” _Structural Health Monitoring_ , vol. 0, no. 0, p. 14759217221075241, 0. [Online]. Available: https://doi.org/10.1177/14759217221075241
  * [41] G. Strang, _Computational Science and Engineering_. Wellesley-Cambridge Press, 2007.
  * [42] G. Golub and C. Van Loan, _Matrix Computations_ , ser. Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, 2013.
  * [43] F. Chung and S.-T. Yau, “Discrete green’s functions,” _Journal of Combinatorial Theory, Series A_ , vol. 91, no. 1-2, pp. 191–214, 2000.
  * [44] N. Srebro and A. Shraibman, “Rank, trace-norm and max-norm,” in _Proc. of the International Conference on Computational Learning Theory_. Springer, 2005, pp. 545–560.
  * [45] Y. Xu and W. Yin, “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,” _SIAM Journal on imaging sciences_ , vol. 6, no. 3, pp. 1758–1789, 2013.
  * [46] C. Malherbe and N. Vayatis, “Global optimization of lipschitz functions,” in _International Conference on Machine Learning_. PMLR, 2017, pp. 2314–2323.
  * [47] A. V. Oppenheim and R. W. Schafer, _Digital signal processing_ , ser. Prentice-Hall international editions. Prentice-Hall, 1975. [Online]. Available: https://www.worldcat.org/oclc/01031348
  * [48] G. Strang, “Wavelets and dilation equations: A brief introduction,” _SIAM Rev._ , vol. 31, no. 4, pp. 614–627, 1989. [Online]. Available: https://doi.org/10.1137/1031128
  * [49] C. Dorn, S. Dasari, Y. Yang, C. Farrar, G. Kenyon, P. Welch, and D. Mascareñas, “Efficient full-field vibration measurements and operational modal analysis using neuromorphic event-based imaging,” _Journal of Engineering Mechanics_ , vol. 144, no. 7, 7 2018.
  * [50] Y. Yang, C. Dorn, T. Mancini, Z. Talken, G. Kenyon, C. Farrar, and D. Mascareñas, “Blind identification of full-field vibration modes from video measurements with phase-based video motion magnification,” _Mechanical Systems and Signal Processing_ , vol. 85, pp. 567–590, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0888327016303272
  * [51] H. Ma, Y.-Z. Lin, and Z. Nie, “Physical interpretation of principal component analysis for structural dynamics through string vibration,” _International Journal of Structural Stability and Dynamics_ , 2019.
  * [52] Z. Lai, I. Alzugaray, M. Chli, and E. Chatzi, “Full-field structural monitoring using event cameras and physics-informed sparse identification,” _Mechanical Systems and Signal Processing_ , vol. 145, p. 106905, 2020.
  * [53] S. Haykin, “Radar signal processing,” _IEEE ASSP Magazine_ , vol. 2, no. 2, pp. 2–18, 1985.
  * [54] J. A. Jensen, _Estimation of blood velocities using ultrasound: a signal processing approach_. Cambridge University Press, 1996.
  * [55] Y. Huang, J. Benesty, and J. Chen, _Acoustic MIMO signal processing_. Springer Science & Business Media, 2006.
  * [56] D. S. Jones, “Introduction to wave phenomena. by a. hirose and k. e. lonngren. wiley, 1985. 355 pp. £53.” _Journal of Fluid Mechanics_ , vol. 162, p. 503–503, 1986.
  * [57] H. Akaike, “Information theory and an extension of the maximum likelihood principle,” in _Selected papers of hirotugu akaike_. Springer, 1998, pp. 199–213.
  * [58] Z. Liao, L. Song, P. Chen, Z. Guan, Z. Fang, and K. Li, “An effective singular value selection and bearing fault signal filtering diagnosis method based on false nearest neighbors and statistical information criteria,” _Sensors (Basel, Switzerland)_ , vol. 18, 2018.
  * [59] P. Comon, “Independent component analysis, a new concept?” _Signal Processing_ , vol. 36, no. 3, pp. 287–314, 1994, higher Order Statistics. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0165168494900299
  * [60] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, “On dynamic mode decomposition: Theory and applications,” _Journal of Computational Dynamics_ , vol. 1, no. 2, pp. 391–421, 2014.
  * [61] Z. Wu and N. E. Huang, “Ensemble empirical mode decomposition: a noise-assisted data analysis method,” _Advances in adaptive data analysis_ , vol. 1, no. 01, pp. 1–41, 2009.
  * [62] D. Kim, M. Park, and H.-S. Oh, “Bidimensional statistical empirical mode decomposition,” _IEEE Signal Processing Letters_ , vol. 19, no. 4, pp. 191–194, 2012.
  * [63] C. Damerval, S. Meignen, and V. Perrier, “A fast algorithm for bidimensional emd,” _IEEE Signal Processing Letters_ , vol. 12, no. 10, pp. 701–704, 2005\.
  * [64] R. J. Allemang, A. W. Phillips, and M. R. Allemang, “Application of principal component analysis methods to experimental structural dynamics,” in _Structural Dynamics, Volume 3_. Springer, 2011, pp. 477–498.
  * [65] G. Rilling and P. Flandrin, “One or two frequencies? the empirical mode decomposition answers,” _IEEE transactions on signal processing_ , vol. 56, no. 1, pp. 85–95, 2007.
  * [66] B. R. Noack, W. Stankiewicz, M. Morzyński, and P. J. Schmid, “Recursive dynamic mode decomposition of transient and post-transient wake flows,” _Journal of Fluid Mechanics_ , vol. 809, pp. 843–872, 2016.

|  Harsha Vardhan Tetali has a PhD in Electrical and Computer Engineering from
University of Florida, Gainesville. Prior to this, he received his B.Tech in
Electronics and Communication Engineering from SVNIT, Surat and M.Tech in
Electrical Engineering from IIT Gandhinagar. His research spans signal
processing, machine learning, and optimization theory. He currently works as a
staff engineer at Marvell Technology.
---|---
|  Joel B. Harley (S’05-M’14) received his B.S. degree in Electrical
Engineering from Tufts University in Medford, MA, USA. He received his M.S.
and Ph.D. degrees in Electrical and Computer Engineering from Carnegie Mellon
University in Pittsburgh, PA, USA in 2011 and 2014, respectively. In 2018, he
joined the University of Florida, where he is currently an associate professor
in the Department of Electrical and Computer Engineering. Previously, he was
an assistant professor in the Department of Electrical and Computer
Engineering at the University of Utah. His research interests include
integrating novel signal processing, machine learning, and data science
methods for the analysis of waves and time-series data. Dr. Harley’s awards
and honors include the 2021 Achenbach Medal from the International Workshop on
Structural Health Monitoring, a 2020 IEEE Ultrasonics, Ferroelectrics, and
Frequency Control Society Star Ambassador Award, a 2020 and 2018 Air Force
Summer Faculty Fellowship, a 2017 Air Force Young Investigator Award, a 2014
Carnegie Mellon A. G. Jordan Award (for academic excellence and exceptional
service to the community). He has published more than 90 technical journal and
conference papers, including four best student papers. He is a member of the
Ultrasonics, Ferroelectrics, and Frequency Control Society, a member of the
IEEE Signal Processing Society, and a member of the Acoustical Society of
America.
---|---
|  Benjamin D. Haeffele is an Associate Research Scientist in the Mathematical
Institute for Data Science and the Center for Imaging Science at Johns Hopkins
University. His research interests involve developing theory and algorithms
for processing high-dimensional data at the intersection of machine learning,
optimization, and computer vision. In addition to basic research in data
science he also works on a variety of applications in medicine, microscopy,
and computational imaging. He received his Ph.D. in Biomedical Engineering at
Johns Hopkins University in 2015 and his B.S. in Electrical Engineering from
the Georgia Institute of Technology in 2006.
---|---

##  VII Supplementary Material

###  VII-A Global Optimal Solution is Achievable

######  Proposition 1.

The optimization problem in (14) is a special case of the problem considered
in [38], which shows that a global minimum is achievable.

###### Proof.

The problem considered in [38] is stated in (15). Comparing it with (14) we
have:

$\displaystyle\ell\left(\mathbf{D}\mathbf{X}^{\top}\right)=\tfrac{1}{2}\|\mathbf{Y}-\mathbf{D}\mathbf{X}^{\top}\|_{F}^{2}$
(37)
$\displaystyle\widebar{\theta}\left(\mathbf{d}_{i},\mathbf{x}_{i}\right)=\tfrac{1}{2}\\!\left(\|\mathbf{x}_{i}\|_{F}^{2}+\min_{k_{i}}\mathbf{d}_{i}^{\top}\mathbf{A}(k_{i})\mathbf{d}_{i}\right)$
(38)

Observe that
$\ell(\widehat{\mathbf{Y}})=\tfrac{1}{2}\|\mathbf{Y}-\widehat{\mathbf{Y}}\|^{2}_{F}$
is convex and differentiable w.r.t. $\widehat{\mathbf{Y}}$. To realize that
the optimization problem in (14) is a special case of the problem considered
in [38], it suffices to check that $\widebar{\theta}(\mathbf{d},\mathbf{x})$
satisfies three conditions:

  1. 1.

$\widebar{\theta}(\alpha\mathbf{d},\alpha\mathbf{x})=\alpha^{2}\theta(\mathbf{d},\mathbf{x}),\
\forall(\mathbf{d},\mathbf{x})$ and $\forall\alpha\geq 0$.

For any $\alpha>0$, $\forall(\mathbf{d},\mathbf{x})$ :

$\displaystyle\widebar{\theta}\left(\alpha\mathbf{d},\alpha\mathbf{x}\right)$
$\displaystyle=\|\alpha\mathbf{d}\|^{2}_{2}+\|\alpha\mathbf{x}\|^{2}_{2}+\gamma\min_{k}\|\alpha\mathbf{L}\mathbf{d}+k^{2}\alpha\mathbf{d}\|^{2}_{F}$
$\displaystyle=\alpha^{2}\|\mathbf{d}\|^{2}_{2}+\alpha^{2}\|\mathbf{x}\|^{2}_{2}+\alpha^{2}\gamma\min_{k}\|\mathbf{L}\mathbf{d}+k^{2}\mathbf{d}\|^{2}_{F}$
$\displaystyle=\alpha^{2}\widebar{\theta}(\mathbf{d},\mathbf{x})$ (39)

where we note that scaling $\mathbf{d}$ by $\alpha>0$ does not change the
optimal value of $k$ in the third term, allowing $\alpha^{2}$ to be moved
outside of the norm.

  2. 2.

$\widebar{\theta}(\mathbf{d},\mathbf{x})\geq 0,\
\forall(\mathbf{d},\mathbf{x})$.

All terms in $\widebar{\theta}(\mathbf{d},\mathbf{x})$ are non-negative, thus,
$\forall(\mathbf{d},\mathbf{x})$, we have
$\widebar{\theta}(\mathbf{d},\mathbf{x})\geq 0$.

  3. 3.

For all sequences $(\mathbf{d}^{(n)},\mathbf{x}^{(n)})$ such that
$\|\mathbf{d}^{(n)}(\mathbf{x}^{(n)})^{\top}\|\rightarrow\infty$ then
$\widebar{\theta}(\mathbf{d}^{(n)},\mathbf{x}^{(n)})\rightarrow\infty$.

Here, note that the following is true for all $(\mathbf{d},\mathbf{x})$:

$\begin{split}\|\mathbf{d}\mathbf{x}^{\top}\|_{F}=\|\mathbf{d}\|_{2}\|\mathbf{x}\|_{2}\leq\tfrac{1}{2}(\|\mathbf{d}\|_{2}^{2}+\|\mathbf{x}\|_{2}^{2})\\\
\leq\tfrac{1}{2}\left(\|\mathbf{d}\|_{2}^{2}+\|\mathbf{x}\|_{2}^{2}+\min_{k}\|\mathbf{L}\mathbf{d}+k^{2}\mathbf{d}\|^{2}_{2}\right)\end{split}$
(40)

As a result we have $\forall(\mathbf{d},\mathbf{x})$ that
$\|\mathbf{d}\mathbf{x}^{\top}\|_{F}\leq\widebar{\theta}(\mathbf{d},\mathbf{x})$,
completing the result.

Hence all of the conditions are satisfied for the global optimal solution to
be achievable. ∎

###  VII-B Solution to the Polar Problem

We restate theorem 1 below.

######  Theorem 1.

For (14), the polar problem in (16) is

$\displaystyle\Omega_{\theta}^{\circ}(\mathbf{Z})=$
$\displaystyle\max_{\mathbf{d},\mathbf{x},k}\,\mathbf{d}^{\top}\mathbf{Z}\mathbf{x}$
$\displaystyle\mathrm{s.t.}~{}\mathbf{d}^{\top}\mathbf{A}(k)\mathbf{d}\leq
1,\|\mathbf{x}\|^{2}_{F}\leq 1,\ 0\leq k\leq 2.$

where $\mathbf{A}(k)$ is as defined in (13). Further, if we define $k^{*}$ as

$k^{*}=\operatorname*{arg\,max}_{k\in[0,2]}\|\mathbf{A}(k)^{-1/2}\mathbf{Z}\|_{2}\;.$
(41)

Then the optimal values of $\mathbf{d},\mathbf{x},k$ are given as
$\mathbf{d}^{*}=\mathbf{A}(k^{*})^{-1/2}\widebar{\mathbf{d}}$,
$\mathbf{x}^{*}=\widebar{\mathbf{x}}$, and $k^{*}$. Where
$\widebar{\mathbf{d}}$ and $\widebar{\mathbf{x}}$ are the left and right
singular vectors, respectively, associated with the largest singular value of
$\mathbf{A}(k^{*})^{-1/2}\mathbf{Z}$.

###### Proof.

The polar problem associated with the objective of the form (15) as given in
[38] is:

$\displaystyle\Omega_{\theta}^{\circ}(\mathbf{Z})=\sup_{\mathbf{d},\mathbf{x}}\mathbf{d}^{\top}\mathbf{Z}\mathbf{x}\
\ \textnormal{s.t.}\ \ \widebar{\theta}(\mathbf{d},\mathbf{x})\leq 1$ (42)

For our particular problem, due to the bilinearity between $\mathbf{d}$ and
$\mathbf{x}$ in the objective the above is equivalent to:

$\displaystyle\begin{split}\Omega_{\theta}^{\circ}(\mathbf{Z})&=\sup_{\mathbf{d},\mathbf{x}}\mathbf{d}^{\top}\mathbf{Z}\mathbf{x}\\\
\ \ \textnormal{s.t.}\ \ \|\mathbf{x}\|^{2}_{2}&\leq
1,\|\mathbf{d}\|^{2}_{2}+\gamma\min_{k}\|\mathbf{L}\mathbf{d}+k^{2}\mathbf{d}\|^{2}_{2}\leq
1\end{split}$ (43)

Note that this is equivalent to moving the minimization w.r.t. $k$ in the
regularization constraint to a maximization over $k$:

$\displaystyle\begin{split}\Omega_{\theta}^{\circ}(\mathbf{Z})&=\sup_{\mathbf{d},\mathbf{x},k}\mathbf{d}^{\top}\mathbf{Z}\mathbf{x}\\\
\ \ \textnormal{s.t.}\ \ \|\mathbf{x}\|^{2}_{2}&\leq
1,\|\mathbf{d}\|^{2}_{2}+\gamma\|\mathbf{L}\mathbf{d}+k^{2}\mathbf{d}\|^{2}_{2}\leq
1\end{split}$ (44)

Next, note that maximizing w.r.t. $\mathbf{d}$ while holding $\mathbf{x}$ and
$k$ fixed is equivalent to solving a problem of the form:

$\displaystyle\max_{\mathbf{d}}\langle\mathbf{d},\mathbf{Z}\mathbf{x}\rangle\
\ \textnormal{s.t.}\ \ \mathbf{d}^{\top}\mathbf{A}(k_{i})\mathbf{d}\leq 1$
(45)

for some positive definite matrix $\mathbf{A}(k_{i})$. If we make the change
of variables $\widebar{\mathbf{d}}=\mathbf{A}(k_{i})^{1/2}\mathbf{d}$, this
then becomes:

$\displaystyle\max_{\widebar{\mathbf{d}}}\ \ $
$\displaystyle\langle\widebar{\mathbf{d}},\mathbf{A}(k_{i})^{-1/2}\mathbf{Z}\mathbf{x}\rangle\
\ \textnormal{s.t.}\ \ \|\widebar{\mathbf{d}}\|_{2}^{2}\leq
1=\|\mathbf{A}(k_{i})^{-1/2}\mathbf{Z}\mathbf{x}\|_{2}$ (46)

where the optimal $\widebar{\mathbf{d}}$ and $\mathbf{d}$ are obtained at

$\displaystyle\widebar{\mathbf{d}}_{opt}$
$\displaystyle=\frac{\mathbf{A}^{-1/2}\mathbf{Z}\mathbf{x}}{\|\mathbf{A}^{-1/2}\mathbf{Z}\mathbf{x}\|_{2}}$
(47) $\displaystyle\mathbf{d}_{opt}$
$\displaystyle=\mathbf{A}^{-1/2}\widebar{\mathbf{d}}_{opt}=\frac{\mathbf{A}^{-1}\mathbf{Z}\mathbf{x}}{\|\mathbf{A}^{-1/2}\mathbf{Z}\mathbf{x}\|_{2}}$
(48)

For our particular problem, if we make the change of variables
$\widebar{k}=k^{2}$ we have that $\mathbf{A}$ is given by:

$\mathbf{A}(\widebar{k})=(1+\gamma\widebar{k}^{2})\mathbf{I}+\gamma\mathbf{L}^{2}+2\widebar{k}\gamma\mathbf{L}$
(49)

where we have used that $\mathbf{L}$ is a symmetric matrix. If we let
$\mathbf{L}=\Gamma\Lambda\Gamma^{\top}$ be an eigen-decomposition of
$\mathbf{L}$ then we can also represent $\mathbf{A}(\widebar{k})$ and
$\mathbf{A}(\widebar{k})^{-1/2}$ as:

$\displaystyle\mathbf{A}(\widebar{k})$
$\displaystyle=\Gamma\left((1+\gamma\widebar{k}^{2})\mathbf{I}+\gamma\Lambda^{2}+2\widebar{k}\gamma\Lambda\right)\Gamma^{\top}$
(50)
$\displaystyle=\Gamma(\mathbf{I}+\gamma(\widebar{k}\mathbf{I}+\Lambda)^{2})\Gamma^{\top}$
(51) $\displaystyle\mathbf{A}(\widebar{k})^{-1/2}$
$\displaystyle=\Gamma(\mathbf{I}+\gamma(\widebar{k}\mathbf{I}+\Lambda)^{2})^{-1/2}\Gamma^{\top}$
(52)

Now we substitute back into the original polar problem:

$\displaystyle\Omega^{\circ}(\mathbf{Z})$
$\displaystyle=\max_{\mathbf{x},\widebar{k}}\|\mathbf{A}(\widebar{k})^{-1/2}\mathbf{Z}\mathbf{x}\|_{2}\
\ \textnormal{subject to}\ \ \|\mathbf{x}\|_{2}^{2}\leq 1$ (53)
$\displaystyle=\max_{\widebar{k}}\|(\mathbf{A}(\widebar{k})^{-1/2}\mathbf{Z}\|_{2}$
(54)

where $\|\cdot\|_{2}$ denotes the spectral norm (maximum singular value).
Similarly, for a given $\widebar{k}$ the optimal $\mathbf{x}$ is given as the
right singular vector of $\mathbf{A}(\widebar{k})^{-1/2}\mathbf{Z}$ associated
with the largest singular value.

As a result, we can solve the polar by performing a line search over
$\widebar{k}$, then once an optimal $\widebar{k}^{*}$ is found we get
$\mathbf{x}^{*}$ as the largest right singular vector of
$\mathbf{A}(\widebar{k}^{*})^{-1/2}\mathbf{Z}$ and the optimal
$\mathbf{d}^{*}$ from (48) (where $\mathbf{d}_{opt}$ will be the largest left
singular vector of $\mathbf{A}(\widebar{k}^{*})^{-1/2}\mathbf{Z}$ multiplied
by $\mathbf{A}(\widebar{k}^{*})^{-1/2}$).

Now, an upper bound for $\widebar{k}$ can be calculated from the fact that the
optimal $\widebar{k}$ is defined using a minimization problem, i.e.

$\displaystyle\widebar{k}^{*}=\operatorname*{arg\,min}_{\widebar{k}}\|\mathbf{L}\mathbf{d}+\widebar{k}\mathbf{d}\|^{2}$
(55)

So we note for any $\mathbf{d}$,

$\displaystyle\widebar{k}^{*}=-\frac{\mathbf{d}^{\top}\mathbf{L}\mathbf{d}}{\|\mathbf{d}\|_{2}^{2}}$
(56)

which is bounded by the smallest eigenvalue of $\mathbf{L}$ (note that
$\mathbf{L}$ is negative (semi)definite). We cite the literature on
eigenvalues of discrete second derivatives [43] to note that all eigenvalues
of $\mathbf{L}$ (irrespective of the boundary conditions) lie in the range
$[\frac{-4}{\Delta l},0]$, since we specifically chose $\Delta l=1$ (without
loss of generality), we have that all eigenvalues of $\mathbf{L}$ lie in the
range$[-4,0]$.

As a result, we need to only consider $\widebar{k}$ in the range $[0,4]$:

$\widebar{k}^{*}=\operatorname*{arg\,max}_{\widebar{k}\in[0,4]}\|\mathbf{A}(\widebar{k})^{-1/2}\mathbf{Z}\|_{2}$
(57)

∎

###  VII-C Optimization Algorithm is Solvable in Polynomial TIme

We first prove a lemma for Lipschitz constants.

######  Lemma 1.

Given a set of constants $\lambda_{i},\ i=1,2,\ldots$ such that $\forall i$,
$\mu_{\lambda}\leq\lambda_{i}\leq 0$, and a constant $\gamma>0$, let the
function $f$ be defined as:

$f(x)=\max_{i}\frac{1}{\sqrt{1+\gamma(x+\lambda_{i})^{2}}}.$ (58)

Then, over the domain $0\leq x\leq\mu_{x}$ $f$ is Lipschitz continuous with
Lipschitz constant $L_{f}$ bounded as follows:

$\begin{split}L_{f}&\leq\begin{cases}\frac{2}{3\sqrt{3}}\sqrt{\gamma},\qquad\gamma\geq\frac{1}{2\max\\{\mu_{\lambda}^{2},\mu_{x}^{2}\\}}\\\
\gamma\max\\{-\mu_{\lambda},\mu_{x}\\}(1+\gamma\max\\{\mu_{\lambda}^{2},\mu_{x}^{2}\\})^{-\tfrac{3}{2}},\quad\mathrm{otherwise}\end{cases}\\\
&\leq\frac{2}{3\sqrt{3}}\sqrt{\gamma}.\end{split}$ (59)

###### Proof.

First, note that for any two Lipschitz continuous functions $\psi_{a}$ and
$\psi_{b}$, with associated Lipschitz constants $L_{a}$ and $L_{b}$,
respectively, one has that the point-wise maximum of the two functions,
$\psi(x)=\max\\{\psi_{a}(x),\psi_{b}(x)\\}$, is also Lipschitz continuous with
Lipschitz constant bounded by $\max\\{L_{a},L_{b}\\}$. This can be easily seen
by the following two inequalities:

$\begin{split}\psi_{a}(x^{\prime})&\leq\psi_{a}(x)+|\psi_{a}(x^{\prime})-\psi_{a}(x)|\\\
&\leq\psi(x)+|\psi_{a}(x^{\prime})-\psi_{a}(x)|\\\
&\leq\psi(x)+L_{a}|x^{\prime}-x|\\\
\psi_{b}(x^{\prime})&\leq\psi_{b}(x)+|\psi_{b}(x^{\prime})-\psi_{b}(x)|\\\
&\leq\psi(x)+|\psi_{b}(x^{\prime})-\psi_{b}(x)|\\\
&\leq\psi(x)+L_{b}|x^{\prime}-x|\end{split}$ (60)

From this we have:

$\begin{split}\psi(x^{\prime})&=\max\\{\psi_{a}(x^{\prime}),\psi_{b}(x^{\prime})\\}\\\
&\leq\max\\{\psi(x)+L_{a}|x^{\prime}-x|,\psi(x)+L_{b}|x^{\prime}-x|\\}\\\
&=\psi(x)+\max\\{L_{a},L_{b}\\}|x^{\prime}-x|\\\
&\implies\psi(x^{\prime})-\psi(x)\leq\max\\{L_{a},L_{b}\\}|x^{\prime}-x|\end{split}$
(61)

which implies the claim from symmetry.

Now, if we define the functions $g$ and $h$ as

$g(x,\lambda)=\frac{1}{\sqrt{1+\gamma(x+\lambda)^{2}}}\ \ \ \ \ \ \
h(b)=\gamma b(1+\gamma b^{2})^{-\tfrac{3}{2}}$ (62)

we have that the Lipschitz constant of $f$, denoted as $L_{f}$, is bounded as:

$\begin{split}L_{f}&\leq\max_{i}\sup_{x\in[0,u_{x}]}\left|\frac{\partial}{\partial
x}g(x,\lambda_{i})\right|\\\
&\leq\sup_{\lambda\in[\mu_{\lambda},0]}\sup_{x\in[0,\mu_{x}]}\left|\frac{\partial}{\partial
x}g(x,\lambda)\right|\\\
&=\sup_{\lambda\in[\mu_{\lambda},0]}\sup_{x\in[0,\mu_{x}]}\left|-\frac{\gamma(x+\lambda)}{(\sqrt{1+\gamma(x+\lambda)^{2}})^{3}}\right|\\\
&=\sup_{b\in[\mu_{\lambda},\mu_{x}]}\left|h(b)\right|\\\
&=\sup_{b\in[0,\max\\{-\mu_{\lambda},\mu_{x}\\}]}h(b)\end{split}$ (63)

Where the first inequality is from the result above about the Lipschitz
constant of the point-wise maximum of two functions and the simple fact that
the Lipschitz constant of a function is bounded by the maximum magnitude of
its gradient, and the final equality is due to the symmetry of $|h(b)|$ about
the origin. Now, finding the critical points of $h(b)$ for non-negative $b$ we
have:

$\begin{split}h^{\prime}(b)&=\gamma(1+\gamma
b^{2})^{-\tfrac{3}{2}}-3\gamma^{2}b^{2}(1+\gamma b^{2})^{-\tfrac{5}{2}}=0\\\
\implies 3\gamma b^{2}&=(1+\gamma b^{2})\implies
b^{*}=\frac{1}{\sqrt{2\gamma}}\end{split}$ (64)

Note that $h(0)=0$, $h(b)>0$ for all $b>0$, and $h^{\prime}(b)<0$ for all
$b>b^{*}$. As a result, $b^{*}$ will be a maximizer of $h(b)$ if it is
feasible, otherwise the maximum will occur at the extreme point
$b=\max\\{-\mu_{\lambda},\mu_{x}\\}$. From this we have the result:

$\begin{split}L_{f}\\\
&\leq\begin{cases}h(b^{*})=\frac{2}{3\sqrt{3}}\sqrt{\gamma}&\gamma\geq\frac{1}{2\max\\{\mu_{\lambda}^{2},\mu_{x}^{2}\\}}\\\
h(\max\\{-\mu_{\lambda},\mu_{x}\\})&\text{otherwise}\end{cases}\\\
&=\begin{cases}\frac{2}{3\sqrt{3}}\sqrt{\gamma}\qquad\qquad\quad\gamma\geq\frac{1}{2\max\\{\mu_{\lambda}^{2},\mu_{x}^{2}\\}}\\\
\gamma\max\\{-\mu_{\lambda},\mu_{x}\\}(1+\gamma\max\\{\mu_{\lambda}^{2},\mu_{x}^{2}\\})^{-\tfrac{3}{2}}\quad\text{otherwise}\end{cases}\\\
&\leq h(b^{*})=\frac{2}{3\sqrt{3}}\sqrt{\gamma}.\end{split}$ (65)

∎

######  Theorem 2.

The line search in equation (57) over $\widebar{k}$ is Lipschitz continuous
with a Lipschitz constant, $L_{\widebar{k}}$, which is bounded by:

$L_{\widebar{k}}\leq\left[\begin{cases}\frac{2}{3\sqrt{3}}\sqrt{\gamma}\|\mathbf{Z}\|_{2}&\gamma\geq\frac{1}{32}\\\
4\gamma(1+16\gamma)^{-\tfrac{3}{2}}\|\mathbf{Z}\|_{2}&\mathrm{otherwise}\end{cases}\right]\leq\frac{2}{3\sqrt{3}}\sqrt{\gamma}\|\mathbf{Z}\|_{2}$
(66)

###### Proof.

To show Lipschitz continuity, we define the function:

$f_{\mathbf{A}}(\widebar{k})=\|\mathbf{A}(\widebar{k})^{-1/2}\|_{2}$ (67)

and then note the following:

$\begin{split}&\left|\|\mathbf{A}(\widebar{k})^{-1/2}\mathbf{Z}\|_{2}-\|\mathbf{A}(\widebar{k}^{\prime})^{-1/2}\mathbf{Z}\|_{2}\right|\\\
\leq&\left\|\left(\mathbf{A}(\widebar{k})^{-1/2}-\mathbf{A}(\widebar{k}^{\prime})^{-1/2}\right)\mathbf{Z}\right\|_{2}\\\
\leq&\left\|\mathbf{A}(\widebar{k})^{-1/2}-\mathbf{A}(\widebar{k}^{\prime})^{-1/2}\right\|_{2}\left\|\mathbf{Z}\right\|_{2}\\\
\leq&L_{\mathbf{A}}|\widebar{k}-\widebar{k}^{\prime}|\|\mathbf{Z}\|_{2}\end{split}$
(68)

where the first inequality is simply the reverse triangle inequality, the
second inequality is due to the spectral norm being submultiplicative, and
$L_{\mathbf{A}}$ denotes the Lipschitz constant of
$f_{\mathbf{A}}(\widebar{k})$. From the form of
$\mathbf{A}(\widebar{k})^{-1/2}$ in (52) note that we have:

$f_{\mathbf{A}}(\widebar{k})\equiv\|\mathbf{A}(\widebar{k})^{-1/2}\|_{2}=\max_{i}\frac{1}{\sqrt{1+\gamma(\widebar{k}+\Lambda_{i,i})^{2}}}$
(69)

so the result is completed by recalling from our discussion above that
$\Lambda_{i,i}\in[-4,0],\forall i$ and applying Lemma 1.

∎

We also note that the above result implies that we can solve the polar by
first performing a (one-dimensional) line search over $k$, and due to the fact
that the largest singular value of a matrix is a Lipschitz continuous
function, this line search can be solved efficiently by a variety of global
optimization algorithms. For example, we give the following corollary for the
simple algorithm given in [46], and similar results can be obtained for other
algorithms.

######  Corollary 1 (Adapted from Cor 13 in [46]).

For the function $f(\bar{k})=\|\mathbf{A}(\bar{k})^{-1/2}\mathbf{Z}\|_{2}$ as
defined in Theorem 1, if we let $\bar{k}_{1},\ldots\bar{k}_{r}$ denote the
iterates of the LIPO algorithm in [46] then we have $\forall\delta\in(0,1)$
with probability at least $1-\delta$,

$\max_{\bar{k}\in[0,4]}f(\bar{k})-\max_{i=1\ldots
r}f(\bar{k}_{i})\leq\tfrac{8}{3\sqrt{3}}\sqrt{\gamma}\|\mathbf{Z}\|_{2}\frac{\ln(1/\delta)}{r}$
(70)

As a result, we have that the error of the linesearch converges with rate
$\mathcal{O}(1/r)$, where $r$ is the number of function evaluations of
$f(\bar{k})$, to a global optimum, and then, given the optimal value of $k$,
the optimal $(\mathbf{d},\mathbf{x})$ vectors can be computed in closed-form
via singular value decomposition. Taken together this allows us to employ the
Meta-Algorithm defined in Algorithm 1 to solve problem (14).

Note that the above results show that the Lipschitz constant of the line
search in the polar problem (and hence the efficiency of solving the polar
problem) is bounded by a constant times a norm at the point at which we
calculate the polar
$\mathbf{Z}=\nabla\frac{1}{2}\|Y-\mathbf{D}\mathbf{X}^{\top}\|_{F}^{2}=\mathbf{D}\mathbf{X}^{\top}-\mathbf{Y}$.
Note that in Algorithm 1 we evaluated the polar when $\mathbf{D}$ and
$\mathbf{X}$ are first-order stationary points. Here we briefly note that it
is straight-forward to guarantee that the norm of the point at which we
evaluate the polar (the residual of the least squared error term,
$\mathbf{D}\mathbf{X}^{\top}-\mathbf{Y}$) is always globally bounded at first-
order stationary points.

To see this, consider a very simplified problem, where we hold the directions
of $\mathbf{D}$ and $\mathbf{X}$ fixed in (14) and we simply optimize over a
scalar multiplication of $\mathbf{X}\rightarrow\alpha\mathbf{X}$. I.e., we
want to solve:

$\min_{\alpha}\tfrac{1}{2}\|\mathbf{Y}-\mathbf{D}(\alpha\mathbf{X})^{\top}\|_{F}^{2}+\tfrac{\lambda}{2}\|(\alpha\mathbf{X})\|_{F}^{2}+\tfrac{\lambda}{2}\sum_{i=1}^{N}\mathbf{d}_{i}^{\top}\mathbf{A}(k_{i})\mathbf{d}_{i}$
(71)

It is easily shown that the optimal value for $\alpha$ in this case (assuming
we are not in the trivial case where $\mathbf{D}$ or $\mathbf{X}$ are all-
zeros) is given as:

$\alpha^{*}=\frac{\langle\mathbf{D}\mathbf{X}^{\top},\mathbf{Y}\rangle}{\|\mathbf{D}\mathbf{X}^{\top}\|_{F}^{2}+\lambda\|\mathbf{X}\|_{F}^{2}}$
(72)

Evaluating the residual for this optimal scaling of $\mathbf{X}$ gives:

$\begin{split}&\tfrac{1}{2}\|\mathbf{Y}-\mathbf{D}(\alpha^{*}\mathbf{X})^{\top}\|_{F}^{2}\\\
=&\tfrac{1}{2}\||\mathbf{Y}\|_{F}^{2}+\frac{(\langle\mathbf{Y},\mathbf{D}\mathbf{X}^{\top}\rangle)^{2}}{\|\mathbf{D}\mathbf{X}^{\top}\|_{F}^{2}+\lambda\|\mathbf{X}\|_{F}^{2}}\left(\frac{\|\mathbf{D}\mathbf{X}^{\top}\|_{F}^{2}}{\|\mathbf{D}\mathbf{X}^{\top}\|_{F}^{2}+\lambda\|\mathbf{X}\|_{F}^{2}}-1\right)\\\
&\leq\tfrac{1}{2}\|\mathbf{Y}\|_{F}^{2}\end{split}$ (73)

where the final inequality can be seen by noting the term in the final
parenthesis is always non-positive and the leading coefficient is non-
negative. Note that any first-order stationary point will be at an optimal
scaling of $\alpha$ since otherwise we can improve the objective by simply
scaling $\mathbf{X}$, which would imply we are not at a first order stationary
point. As a result, we have that the residual term when we evaluate it in
Algorithm 1 will always be bounded at least by
$\tfrac{1}{2}\|\mathbf{Y}\|_{F}^{2}$, which also trivially implies a bound on
the spectral norm of the residual since we always have
$\|\mathbf{Z}\|_{2}\leq\|\mathbf{Z}\|_{F}$, $\forall\mathbf{Z}$.

###  VII-D Proof of Optimal Step Size

###### Proof.

Let $f(\tau)$ represent the objective function in (24) with everything except
for $\tau$ held fixed:

$\begin{split}f(\tau)&=\tfrac{1}{2}\|\mathbf{Y}-\mathbf{D}_{\tau}\mathbf{X}_{\tau}^{\top}\|_{F}^{2}+\tfrac{\lambda}{2}\sum_{i=1}^{N+1}(\mathbf{d}_{\tau})_{i}^{\top}\mathbf{A}_{\tau}(k_{i})(\mathbf{d}_{\tau})_{i}\end{split}$
(74)

Observe that from the solution to the Polar problem we have by construction
that the regularization term $\theta(\mathbf{d}^{*},\mathbf{x}^{*})=1$, so
combined with the positive homogeneity of $\theta$ we have have that
minimizing $f(\tau)$ w.r.t. $\tau$ is equivalent to solving:

$\min_{\tau\geq
0}\tfrac{1}{2}\|\mathbf{Y}-\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top}-\tau^{2}\mathbf{d}^{*}(\mathbf{x}^{*})^{\top}\|_{F}^{2}+\lambda\tau^{2}$
(75)

Taking the gradient w.r.t. $\tau^{2}$ and solving for 0 gives:

$(\tau^{*})^{2}=\frac{\langle\mathbf{Y}-\widetilde{\mathbf{D}}\widetilde{\mathbf{X}}^{\top},\mathbf{d}^{*}(\mathbf{x}^{*})^{\top}\rangle-\lambda}{\|\mathbf{d}^{*}(\mathbf{x}^{*})^{\top}\|_{F}^{2}}$
(76)

The result is completed by noting that the numerator is guaranteed to be
strictly positive due to the fact that the Polar solution has value strictly
greater that 1. ∎