2}_{1}\lambda^{2}_{1}+\zeta^{\langle I\rangle
2}_{2}\lambda^{2}_{2}}{\zeta^{\langle I\rangle
2}_{1}\lambda^{3}_{1}+\zeta^{\langle I\rangle
2}_{2}\lambda^{3}_{2}}\lambda_{b})\zeta^{\langle
I\rangle}_{b}\boldsymbol{v}_{b}\rVert_{2}}{\lVert\zeta^{\langle
I\rangle}_{1}\boldsymbol{v}_{1}+\zeta^{\langle
I\rangle}_{2}\boldsymbol{v}_{2}\rVert_{2}}$
$\displaystyle=\frac{\lVert\sum_{b}(1-\frac{\iota^{2}+\tau^{\langle I\rangle
2}}{\lambda_{2}(\iota^{3}+\tau^{\langle I\rangle
2})}\lambda_{b})\zeta^{\langle
I\rangle}_{b}\boldsymbol{v}_{b}\rVert_{2}}{\lVert\zeta^{\langle
I\rangle}_{1}\boldsymbol{v}_{1}+\zeta^{\langle
I\rangle}_{2}\boldsymbol{v}_{2}\rVert_{2}}$
$\displaystyle=\frac{\lVert(1-\frac{\iota^{2}+\tau^{\langle I\rangle
2}}{\lambda_{2}(\iota^{3}+\tau^{\langle I\rangle
2})}\lambda_{1})\zeta^{\langle
I\rangle}_{1}\boldsymbol{v}_{1}+(1-\frac{\iota^{2}+\tau^{\langle I\rangle
2}}{\lambda_{2}(\iota^{3}+\tau^{\langle I\rangle
2})}\lambda_{2})\zeta^{\langle
I\rangle}_{2}\boldsymbol{v}_{2}\rVert_{2}}{\lVert\zeta^{\langle
I\rangle}_{1}\boldsymbol{v}_{1}+\zeta^{\langle
I\rangle}_{2}\boldsymbol{v}_{2}\rVert_{2}}$
$\displaystyle=\frac{\lVert(\zeta^{\langle
I\rangle}_{1}\boldsymbol{v}_{1}+\zeta^{\langle
I\rangle}_{2}\boldsymbol{v}_{2})-(\iota\frac{\iota^{2}+\tau^{\langle I\rangle
2}}{(\iota^{3}+\tau^{\langle I\rangle 2})})\zeta^{\langle
I\rangle}_{1}\boldsymbol{v}_{1}-(\frac{\iota^{2}+\tau^{\langle I\rangle
2}}{(\iota^{3}+\tau^{\langle I\rangle 2})})\zeta^{\langle
I\rangle}_{2}\boldsymbol{v}_{2}\rVert_{2}}{\lVert\zeta^{\langle
I\rangle}_{1}\boldsymbol{v}_{1}+\zeta_{i2}\boldsymbol{v}_{2}\rVert_{2}}$
$\displaystyle=\frac{\lVert(\zeta^{\langle
I\rangle}_{1}\boldsymbol{v}_{1}+\zeta^{\langle
I\rangle}_{2}\boldsymbol{v}_{2})-\frac{\iota^{2}+\tau^{\langle I\rangle
2}}{(\iota^{3}+\tau^{\langle I\rangle 2})}(\iota\zeta^{\langle
I\rangle}_{1}\boldsymbol{v}_{1}+\zeta^{\langle
I\rangle}_{2}\boldsymbol{v}_{2})\rVert_{2}}{\lVert\zeta^{\langle
I\rangle}_{1}\boldsymbol{v}_{1}+\zeta^{\langle
I\rangle}_{2}\boldsymbol{v}_{2}\rVert_{2}}$
$\displaystyle=\frac{\lVert(\boldsymbol{v}_{1}+\tau^{\langle
I\rangle}\boldsymbol{v}_{2})-\frac{\iota^{2}+\tau^{\langle I\rangle
2}}{(\iota^{3}+\tau^{\langle I\rangle
2})}(\iota\boldsymbol{v}_{1}+\tau^{\langle
I\rangle}\boldsymbol{v}_{2})\rVert_{2}}{\lVert\boldsymbol{v}_{1}+\tau^{\langle
I\rangle}\boldsymbol{v}_{2}\rVert_{2}}.$

The worst convergence (i.e. an upper bound) is achieved when
$\iota=\tau^{\langle I\rangle}$. Hence, we write the upper bound on
$\omega^{\langle I\rangle}$ as follows:

$\displaystyle\omega^{\langle I\rangle}$
$\displaystyle=\frac{\lVert(\boldsymbol{v}_{1}+\tau^{\langle
I\rangle}\boldsymbol{v}_{2})-\frac{\iota^{2}+\tau^{\langle I\rangle
2}}{(\iota^{3}+\tau^{\langle I\rangle
2})}(\iota\boldsymbol{v}_{1}+\tau^{\langle
I\rangle}\boldsymbol{v}_{2})\rVert_{2}}{\lVert\boldsymbol{v}_{1}+\tau^{\langle
I\rangle}\boldsymbol{v}_{2}\rVert_{2}}$ (94)
$\displaystyle\leq\frac{\lVert(\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2})-\frac{\iota^{2}+\iota^{2}}{(\iota^{3}+\iota^{2})}(\iota\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2})\rVert_{2}}{\lVert\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2}\rVert_{2}}$
$\displaystyle=\frac{\lVert(\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2})-\frac{2\iota^{2}}{\iota^{2}(1+\iota)}(\iota\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2})\rVert_{2}}{\lVert\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2}\rVert_{2}}$
$\displaystyle=\frac{\lVert(\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2})-\frac{2\iota}{1+\iota}(\boldsymbol{v}_{1}+\boldsymbol{v}_{2})\rVert_{2}}{\lVert\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2}\rVert_{2}}$
$\displaystyle=\frac{\lVert\frac{1-\iota}{\iota+1}\boldsymbol{v}_{1}+\frac{\iota(\iota-1)}{\iota+1}\boldsymbol{v}_{2}\rVert_{2}}{\lVert\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2}\rVert_{2}}$
$\displaystyle=\frac{\frac{\iota-1}{\iota+1}\lVert\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2}\rVert_{2}}{\lVert\boldsymbol{v}_{1}+\iota\boldsymbol{v}_{2}\rVert_{2}}$
$\displaystyle=\frac{\iota-1}{\iota+1}.$

Hence,

$\lVert\boldsymbol{\beta}^{\langle
I+1\rangle}-\boldsymbol{\beta}^{\\#}\rVert_{2}\leq\left(\frac{\iota-1}{\iota+1}\right)\lVert\boldsymbol{\beta}^{\langle
I\rangle}-\boldsymbol{\beta}^{\\#}\rVert_{2}.$ (95)

We denote $\rho_{\text{SD}}\coloneqq\frac{\iota-1}{\iota+1}$ (which is a
faster rate than $\rho_{\text{GD}}$), and write

$\lVert\boldsymbol{\beta}^{\langle
I\rangle}-\boldsymbol{\beta}^{\\#}\rVert_{2}=\mathcal{O}(\rho_{\text{SD}}^{I}).$
(96)

This bound also holds for $D>2$; the proof steps are similar (e.g. see Section
9.2 of Shewchuk et al. (1994)).

####  E.4.2 Jacobian Convergence

We now analyze the convergence rate of the Jacobian. Convergence studies in
this context date back to the seminal work of gilbert1992automatic. We drop
the subscript $\boldsymbol{\theta}$ when referring to $(j,k)$-entry of
$\mathbb{A}_{\boldsymbol{\theta}}$ for ease of notation. Given matrix
$\mathbb{A}$, we have

$\partial\lambda_{j}=\boldsymbol{v}_{j}^{\top}\partial\mathbb{A}_{\boldsymbol{\theta}}\boldsymbol{v}_{j}$
(97)
$\partial\boldsymbol{v}_{j}=(\lambda_{j}\mathbb{I}-\mathbb{A}_{\boldsymbol{\theta}})^{{\dagger}}\partial\mathbb{A}_{\boldsymbol{\theta}}\boldsymbol{v}_{j}$
(98)
$\frac{\partial\mathbb{A}_{\boldsymbol{\theta}}}{\partial\mathbb{A}_{jk}}=\mathbb{H}^{jk},$
(99)

where ${\dagger}$ denotes pseudo inverse and $\mathbb{H}$ is a zero matrix
except at the $(j,k)$-entry, which is $1$. The above holds when the
eigenvalues and vectors are distinct. We denote the Jacobian error by
$\mathbb{B}^{\langle I\rangle}:=\boldsymbol{J}^{\langle
I\rangle}-\boldsymbol{J}^{\\#}=\frac{\partial\boldsymbol{\beta}^{\langle
I\rangle}}{\partial\boldsymbol{\theta}}-\frac{\partial\boldsymbol{\beta}^{\\#}}{\partial\boldsymbol{\theta}}$.
Given the expression $\boldsymbol{\beta}^{\langle
I\rangle}-\boldsymbol{\beta}^{\\#}=\sum_{b}\zeta^{\langle
I\rangle}_{b}\boldsymbol{v}_{b}$, we first focus on writing the Jacobian error
with respect to $\mathbb{A}_{jk}$, denoted by
$\mathbb{B}_{\mathbb{A}_{jk}}^{\langle
I\rangle}=\frac{\partial\boldsymbol{e}^{\langle
I\rangle}}{\partial\mathbb{A}_{jk}}$. Then, we use this Jacobian error to find
the bound on the Jacobian with respect to $\boldsymbol{\theta}$ as
$\mathbb{B}^{\langle
I\rangle}=\sum_{jk}\frac{\partial\mathbb{A}_{jk}}{\partial\boldsymbol{\theta}}\frac{\partial\boldsymbol{e}^{\langle
I\rangle}}{\partial\mathbb{A}_{jk}}$.

##### Gradient Descent

For gradient descent, we start with the recursion below,

$\displaystyle\boldsymbol{\beta}^{\langle
I+1\rangle}-\boldsymbol{\beta}^{\\#}$
$\displaystyle=(\mathbb{I}-\alpha\mathbb{A}_{\boldsymbol{\theta}})(\boldsymbol{\beta}^{\langle
I\rangle}-\boldsymbol{\beta}^{\\#}).$ (100)

Then, we take the derivative with respect to $\mathbb{A}_{jk}$,

$\displaystyle\mathbb{B}^{\langle I+1\rangle}_{\mathbb{A}_{jk}}$
$\displaystyle=(\mathbb{I}-\alpha\mathbb{A}_{\boldsymbol{\theta}})\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}+\frac{\partial(\mathbb{I}-\alpha\mathbb{A}_{\boldsymbol{\theta}})}{\partial\mathbb{A}_{jk}}\boldsymbol{e}^{\langle
I\rangle}$ (101)
$\displaystyle=(\mathbb{I}-\alpha\mathbb{A}_{\boldsymbol{\theta}})\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}-\alpha\mathbb{H}^{jk}\boldsymbol{e}^{\langle
I\rangle}.$

We use the eigendecomposition of
$\mathbb{A}_{\boldsymbol{\theta}}=\mathbb{V}\mathbb{\Sigma}\mathbb{V}^{\top}$
(i.e. $\mathbb{V}$ is the eigenvectors of $\mathbb{A}_{\boldsymbol{\theta}}$,
and $\mathbb{\Sigma}$ is a diagonal matrix of $\mathbb{A}$ eigenvalues.This
leads to the bound

$\displaystyle\lVert\mathbb{B}^{\langle
I+1\rangle}_{\mathbb{A}_{jk}}\rVert_{2}$
$\displaystyle\leq\lVert\mathbb{I}-\alpha\mathbb{A}_{\boldsymbol{\theta}}\rVert_{2}\lVert\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}\rVert_{2}+\lVert\alpha\mathbb{H}^{jk}\boldsymbol{e}^{\langle
I\rangle}\rVert_{2}$ (102)
$\displaystyle\leq\rho_{\text{GD}}\lVert\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}\rVert_{2}+\lVert\alpha\mathbb{H}^{jk}\rVert_{2}\lVert\boldsymbol{e}^{\langle
I\rangle}\rVert_{2}.$

Unrolling the recursion gives us

$\displaystyle\lVert\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}\rVert_{2}=\mathcal{O}(I\rho_{\text{GD}}^{I}).$
(103)

Given $\mathbb{B}^{\langle
I\rangle}=\sum_{jk}\frac{\partial\mathbb{A}_{jk}}{\partial\boldsymbol{\theta}}\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}$, for gradient descent, we have

$\displaystyle\lVert\mathbb{B}^{\langle
I\rangle}\rVert_{2}=\mathcal{O}(I\rho_{\text{GD}}^{I}).$ (104)

##### Steepest Descent

For steepest descent, we analyze Jacobian convergence under two different
scenarios: (a) the gradient with respect to $\boldsymbol{\theta}$ is _not_
propagated through the adaptive step size $\alpha^{\langle I\rangle}$, and (b)
the gradient is propagated through $\alpha^{\langle I\rangle}$ (i.e. a more
sophisticated scenario). In both of these scenarios, we show that the Jacobian
converges at the same asymptotic rate. We express the Jacobian error as

$\displaystyle\mathbb{B}_{\mathbb{A}_{jk}}^{\langle I\rangle}$
$\displaystyle=\sum_{b}\zeta^{\langle
I\rangle}_{b}\frac{\partial\boldsymbol{v}_{b}}{\partial\mathbb{A}_{jk}}$ (105)
$\displaystyle=\sum_{b}\zeta^{\langle
I\rangle}_{b}(\lambda_{b}\boldsymbol{I}-\mathbb{A}_{\boldsymbol{\theta}})^{{\dagger}}\mathbb{H}^{jk}\boldsymbol{v}_{b}$
$\displaystyle=\sum_{b}\zeta^{\langle
I\rangle}_{b}(\lambda_{b}\boldsymbol{I}-\mathbb{V}\mathbb{\Sigma}\mathbb{V}^{\top})^{{\dagger}}\mathbb{H}^{jk}\boldsymbol{v}_{b}$
$\displaystyle=\sum_{b}\zeta^{\langle
I\rangle}_{b}(\mathbb{V}\mathbb{\Lambda}_{b}\mathbb{V}^{\top})\mathbb{H}^{jk}\boldsymbol{v}_{b}$
$\displaystyle=\sum_{b}\sum_{d}\zeta^{\langle
I\rangle}_{b}(q_{bd}\boldsymbol{v}_{jd}\boldsymbol{v}_{bk})\boldsymbol{v}_{d},$

where $\mathbb{\Lambda}_{b}$ is a diagonal matrix with
$q_{bd}:=\frac{1}{\lambda_{b}-\lambda_{d}}$ is its $d$-th diagonal entry, and
$0$ on $b$-th diagonal entry. We substitute this expression into the Jacobian
recursion,

$\displaystyle\mathbb{B}^{\langle I+1\rangle}_{\partial\mathbb{A}_{jk}}$
$\displaystyle=(\mathbb{I}-\alpha^{\langle
I\rangle}\mathbb{A}_{\boldsymbol{\theta}})\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}+\frac{\partial(\mathbb{I}-\alpha^{\langle
I\rangle}\mathbb{A}_{\boldsymbol{\theta}})}{\partial\mathbb{A}_{jk}}\boldsymbol{e}^{\langle
I\rangle}.$ (106)

_Scenario (a)._

$\displaystyle\mathbb{B}^{\langle I+1\rangle}_{\mathbb{A}_{jk}}$
$\displaystyle=(\mathbb{I}-\alpha^{\langle
I\rangle}\mathbb{A}_{\boldsymbol{\theta}})\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}$ (107)
$\displaystyle=(\mathbb{I}-\alpha^{\langle
I\rangle}\mathbb{A}_{\boldsymbol{\theta}})\sum_{b}\sum_{d}\zeta^{\langle
I\rangle}_{b}(q_{bd}\boldsymbol{v}_{jd}\boldsymbol{v}_{bk})\boldsymbol{v}_{d}-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}$
$\displaystyle=\sum_{b}\sum_{d}\zeta^{\langle
I\rangle}_{b}(q_{bd}\boldsymbol{v}_{jd}\boldsymbol{v}_{bk})(\mathbb{I}-\alpha^{\langle
I\rangle}\mathbb{A}_{\boldsymbol{\theta}})\boldsymbol{v}_{d}-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}$
$\displaystyle=\sum_{b}\sum_{d}\zeta^{\langle
I\rangle}_{b}(q_{bd}\boldsymbol{v}_{jd}\boldsymbol{v}_{bk})(1-\alpha^{\langle
I\rangle}\lambda_{d})\boldsymbol{v}_{d}-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}$
$\displaystyle=\sum_{b}\sum_{d}\zeta^{\langle
I\rangle}_{b}(q_{bd}\boldsymbol{v}_{jd}\boldsymbol{v}_{bk})(1-\frac{\sum_{o}\zeta^{\langle
I\rangle 2}_{o}\lambda^{2}_{o}}{\sum_{o}\zeta^{\langle I\rangle
2}_{o}\lambda^{3}_{o}}\lambda_{d})\boldsymbol{v}_{d}-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}.$

Similar to the analysis for the solver error, we consider the case when $D=2$.
We have

$\displaystyle\mathbb{B}^{\langle I+1\rangle}_{\mathbb{A}_{jk}}$
$\displaystyle=\sum_{b}\sum_{d}\zeta^{\langle
I\rangle}_{b}(q_{bd}\boldsymbol{v}_{jd}\boldsymbol{v}_{bk})\boldsymbol{v}_{d}(1-\frac{\sum_{o}\zeta^{\langle
I\rangle 2}_{o}\lambda^{2}_{o}}{\sum_{o}\zeta^{\langle I\rangle
2}_{o}\lambda^{3}_{o}}\lambda_{d})-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}$ (108)
$\displaystyle=\sum_{b}\sum_{d}\zeta^{\langle
I\rangle}_{b}(q_{bd}\boldsymbol{v}_{jd}\boldsymbol{v}_{bk})\boldsymbol{v}_{d}(1-\frac{\sum_{o}\zeta^{\langle
I\rangle 2}_{o}\lambda^{2}_{o}}{\sum_{o}\zeta^{\langle I\rangle
2}_{o}\lambda^{3}_{o}}\lambda_{d})-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}$
$\displaystyle=\sum_{b}\left(\zeta^{\langle
I\rangle}_{b}(q_{b1}\boldsymbol{v}_{j1}\boldsymbol{v}_{bk})\boldsymbol{v}_{1}(1-\frac{\zeta^{\langle
I\rangle 2}_{1}\lambda^{2}_{1}+\zeta^{\langle I\rangle
2}_{2}\lambda^{2}_{2}}{\zeta^{\langle I\rangle
2}_{1}\lambda^{3}_{1}+\zeta^{\langle I\rangle
2}_{2}\lambda^{3}_{2}}\lambda_{1})+\zeta^{\langle
I\rangle}_{b}(q_{b2}\boldsymbol{v}_{j2}\boldsymbol{v}_{bk})\boldsymbol{v}_{2}(1-\frac{\zeta^{\langle
I\rangle 2}_{1}\lambda^{2}_{1}+\zeta^{\langle I\rangle
2}_{2}\lambda^{2}_{2}}{\zeta^{\langle I\rangle
2}_{1}\lambda^{3}_{1}+\zeta^{\langle I\rangle
2}_{2}\lambda^{3}_{2}}\lambda_{2})\right)-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}$
$\displaystyle=\sum_{b}\left(\zeta^{\langle
I\rangle}_{b}(q_{b1}\boldsymbol{v}_{j1}\boldsymbol{v}_{bk})\boldsymbol{v}_{1}(1-\frac{\iota^{2}+\tau^{\langle
I\rangle 2}}{\lambda_{2}(\iota^{3}+\tau^{\langle I\rangle
2})}\lambda_{1})+\zeta^{\langle
I\rangle}_{b}(q_{b2}\boldsymbol{v}_{j2}\boldsymbol{v}_{bk})\boldsymbol{v}_{2}(1-\frac{\iota^{2}+\tau^{\langle
I\rangle 2}}{\lambda_{2}(\iota^{3}+\tau^{\langle I\rangle
2})}\lambda_{2})\right)-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}$
$\displaystyle=\sum_{b}\left(\zeta^{\langle
I\rangle}_{b}(q_{b1}\boldsymbol{v}_{j1}\boldsymbol{v}_{bk})\boldsymbol{v}_{1}(1-\frac{\iota^{2}+\tau^{\langle
I\rangle 2}}{\iota^{3}+\tau^{\langle I\rangle 2}}\iota)+\zeta^{\langle
I\rangle}_{b}(q_{b2}\boldsymbol{v}_{j2}\boldsymbol{v}_{bk})\boldsymbol{v}_{2}(1-\frac{\iota^{2}+\tau^{\langle
I\rangle 2}}{\iota^{3}+\tau^{\langle I\rangle 2}})\right)-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}.$

Now, for the upper bound, we set $\iota=\tau^{\langle I\rangle}$

$\displaystyle\lVert\mathbb{B}^{\langle
I+1\rangle}_{\mathbb{A}_{jk}}\rVert_{2}$
$\displaystyle=\lVert\sum_{b}\left(\zeta^{\langle
I\rangle}_{b}(q_{b1}\boldsymbol{v}_{j1}\boldsymbol{v}_{bk})\boldsymbol{v}_{1}(1-\frac{\iota^{2}+\tau^{\langle
I\rangle 2}}{\iota^{3}+\tau^{\langle I\rangle 2}}\iota)+\zeta^{\langle
I\rangle}_{b}(q_{b2}\boldsymbol{v}_{j2}\boldsymbol{v}_{bk})\boldsymbol{v}_{2}(1-\frac{\iota^{2}+\tau^{\langle
I\rangle 2}}{\iota^{3}+\tau^{\langle I\rangle 2}})\right)-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}\rVert_{2}$ (109)
$\displaystyle\leq\lVert\sum_{b}\left(\zeta^{\langle
I\rangle}_{b}(q_{b1}\boldsymbol{v}_{j1}\boldsymbol{v}_{bk})\boldsymbol{v}_{1}(1-\frac{\iota^{2}+\tau^{\langle
I\rangle 2}}{\iota^{3}+\tau^{\langle I\rangle 2}}\iota)+\zeta^{\langle
I\rangle}_{b}(q_{b2}\boldsymbol{v}_{j2}\boldsymbol{v}_{bk})\boldsymbol{v}_{2}(1-\frac{\iota^{2}+\tau^{\langle
I\rangle 2}}{\iota^{3}+\tau_{i}^{2}})\right)\rVert_{2}+\lVert\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}\rVert_{2}$
$\displaystyle\leq\lVert\sum_{b}\left(\zeta^{\langle
I\rangle}_{b}(q_{b1}\boldsymbol{v}_{j1}\boldsymbol{v}_{bk})\boldsymbol{v}_{1}(1-\frac{2}{1+\iota}\iota)+\zeta^{\langle
I\rangle}_{b}(q_{b2}\boldsymbol{v}_{j2}\boldsymbol{v}_{bk})\boldsymbol{v}_{2}(1-\frac{2}{1+\iota})\right)\rVert_{2}+\lVert\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}\rVert_{2}$
$\displaystyle\leq\lVert\sum_{b}\left(\zeta^{\langle
I\rangle}_{b}(q_{b1}\boldsymbol{v}_{j1}\boldsymbol{v}_{bk})\boldsymbol{v}_{1}(\frac{1-\iota}{\iota+1})+\zeta^{\langle
I\rangle}_{b}(q_{b2}\boldsymbol{v}_{j2}\boldsymbol{v}_{bk})\boldsymbol{v}_{2}(\frac{\iota-1}{\iota+1})\right)\rVert_{2}+\lVert\alpha_{i}\mathbb{H}^{jk}\boldsymbol{e}^{\langle
I\rangle}\rVert_{2}$
$\displaystyle=\lVert\frac{\iota-1}{\iota+1}\sum_{b}\left(\zeta^{\langle
I\rangle}_{b}(q_{b1}\boldsymbol{v}_{j1}\boldsymbol{v}_{bk})\boldsymbol{v}_{1}+\zeta^{\langle
I\rangle}_{b}(q_{b2}\boldsymbol{v}_{j2}\boldsymbol{v}_{bk})\boldsymbol{v}_{2}\right)\rVert_{2}+\lVert\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}\rVert_{2}$
$\displaystyle=\frac{\iota-1}{\iota+1}\lVert\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}\rVert_{2}+\lVert\alpha^{\langle
I\rangle}\mathbb{H}^{jk}\boldsymbol{e}^{\langle I\rangle}\rVert_{2}.$

Unrolling the recursion, we have

$\lVert\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}\rVert=\mathcal{O}(I\rho_{\text{SD}}^{I}).$ (110)

Hence,

$\lVert\mathbb{B}^{\langle
I\rangle}\rVert_{2}=\mathcal{O}(I\rho_{\text{SD}}^{I}).$ (111)

_Scenario (b)._ In scenario (a), we have $\frac{(\mathbb{I}-\alpha^{\langle
I\rangle}\mathbb{A}_{\boldsymbol{\theta}})}{\partial\mathbb{A}_{jk}}=-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}$ with bounded norm of $\alpha^{\langle I\rangle}$. In
scenario (b), we have

$\frac{\partial(\mathbb{I}-\alpha^{\langle
I\rangle}\mathbb{A}_{\boldsymbol{\theta}})}{\partial\mathbb{A}_{jk}}=\frac{\partial\mathbb{A}_{\boldsymbol{\theta}}}{\partial\mathbb{A}_{jk}}\frac{\partial(\mathbb{I}-\alpha^{\langle
I\rangle}\mathbb{A}_{\boldsymbol{\theta}})}{\partial\mathbb{A}_{\boldsymbol{\theta}}}+\frac{\partial\alpha^{\langle
I\rangle}}{\partial\mathbb{A}_{jk}}\frac{\partial(\mathbb{I}-\alpha^{\langle
I\rangle}\mathbb{A}_{\boldsymbol{\theta}})}{\partial\alpha^{\langle
I\rangle}}=-\alpha^{\langle
I\rangle}\mathbb{H}^{jk}-\mathbb{A}_{\boldsymbol{\theta}}\frac{\partial\alpha^{\langle
I\rangle}}{\partial\mathbb{A}_{jk}}.$ (112)

This leads to the following,

$\displaystyle\mathbb{B}^{\langle I+1\rangle}_{\mathbb{A}_{jk}}$
$\displaystyle=(\mathbb{I}-\alpha^{\langle
I\rangle}\mathbb{A}_{\boldsymbol{\theta}})\mathbb{B}^{\langle
I\rangle}_{\mathbb{A}_{jk}}-(\alpha^{\langle
I\rangle}\mathbb{H}^{jk}+\mathbb{A}\frac{\partial\alpha^{\langle
I\rangle}}{\partial\mathbb{A}_{jk}})\boldsymbol{e}^{\langle I\rangle}.$ (113)

It remains to be shown that the norm of
$\mathbb{A}\frac{\partial\alpha^{\langle I\rangle}}{\partial\mathbb{A}_{jk}}$
can be upper-bounded. Then, the Jacobian error is upper-bonded by the same
order of convergence as in scenario (a),

$\lVert\mathbb{A}\frac{\partial\alpha^{\langle
I\rangle}}{\partial\mathbb{A}_{jk}}\rVert_{2}\leq\lVert\mathbb{A}_{\boldsymbol{\theta}}\rVert_{2}\lVert\frac{\partial\alpha^{\langle
I\rangle}}{\partial\mathbb{A}_{jk}}\rVert_{2}=\lambda_{1}\lVert\frac{\partial\alpha^{\langle
I\rangle}}{\partial\mathbb{A}_{jk}}\rVert_{2}.$ (114)

We write $\alpha^{\langle I\rangle}$ in terms of eigenvalues of
$\mathbb{A}_{\boldsymbol{\theta}}$ (i.e. $\alpha^{\langle
I\rangle}=\frac{\sum_{o}\zeta^{\langle I\rangle
2}_{o}\lambda^{2}_{o}}{\sum_{b}\zeta^{\langle I\rangle
2}_{b}\lambda^{3}_{b}}$) and take the derivative, i.e.

$\frac{\partial\alpha^{\langle
I\rangle}}{\partial\mathbb{A}_{jk}}=\frac{\left(\sum_{o}2\zeta^{\langle
I\rangle
2}_{o}\lambda_{o}\boldsymbol{v}_{o}^{\top}\mathbb{H}^{jk}\boldsymbol{v}_{o}\right)\left(\sum_{b}\zeta^{\langle
I\rangle 2}_{b}\lambda^{3}_{b}\right)-\left(\sum_{o}\zeta^{\langle I\rangle
2}_{o}\lambda^{2}_{o}\right)\left(\sum_{b}3\zeta^{\langle I\rangle
2}_{b}\lambda^{2}_{b}\boldsymbol{v}_{b}^{\top}\mathbb{H}^{jk}\boldsymbol{v}_{b}\right)}{(\sum_{m}\zeta^{\langle
I\rangle 2}_{m}\lambda^{3}_{m})^{2}}.$ (115)

Given the above, there exist a constant that bounds the norm of this
derivative. We denote this constant by $\lVert\frac{\partial\alpha^{\langle
I\rangle}}{\partial\mathbb{A}_{jk}}\rVert_{2}\leq M_{\alpha_{\text{div}}}$.

##  Appendix F Additional Experimental Details

###  F.1 Parameter Recovery for Noisy AR Models

####  F.1.1 Form of Prior Precision $\mathbb{\Gamma}_{\boldsymbol{\theta}}$

From (26), let us define $\boldsymbol{z}_{\leq P}:=\\{z_{1},\ldots,z_{P}\\}$,
$\boldsymbol{z}_{>P}:=\\{z_{P+1},\ldots,z_{D}\\}$ and
$\boldsymbol{w}:=\\{w_{P+1},\ldots,w_{D}\\}$. Then,

$\displaystyle\begin{bmatrix}\mathbb{I}&\mathbb{0}\\\
\mathbb{H}&\mathbb{L}\end{bmatrix}\begin{bmatrix}\boldsymbol{z}_{\leq P}\\\
\boldsymbol{z}_{>P}\end{bmatrix}=\begin{bmatrix}\boldsymbol{z}_{\leq P}\\\
\boldsymbol{w}\end{bmatrix},$ (116)

where $\mathbb{H}\in\mathbb{R}^{(D-P)\times P}$ and
$\mathbb{L}\in\mathbb{R}^{(D-P)\times(D-P)}$ such that

$\displaystyle\mathbb{H}:=\begin{bmatrix}-\phi_{P}&-\phi_{P-1}&\ldots&-\phi_{1}\\\
0&-\phi_{P}&\ldots&-\phi_{2}\\\ 0&0&\ldots&\vdots\\\ 0&0&\ldots&-\phi_{P}\\\
0&0&\ldots&0\\\ \vdots&\vdots&\vdots&\vdots\\\ 0&0&\ldots&0\end{bmatrix}$
$\displaystyle\mathbb{L}:=\begin{bmatrix}1&0&0&\ldots&0\\\
-\phi_{1}&1&0&\ldots&0\\\ \vdots&\vdots&\vdots&\ldots&\vdots\\\
-\phi_{P-1}&-\phi_{P-2}&-\phi_{P-3}&\ldots&0\\\
-\phi_{P}&-\phi_{P-1}&-\phi_{P-2}&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\ 0&0&0&\ldots&1\end{bmatrix}.$ (117)

Observe that $\begin{bmatrix}\boldsymbol{z}_{\leq P}\\\
\boldsymbol{z}_{>P}\end{bmatrix}$ is a multivariate Gaussian with mean
$\boldsymbol{0}$ and inverse-covariance
$\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ (i.e. our object of interest).
Similarly, $\begin{bmatrix}\boldsymbol{z}_{\leq P}\\\
\boldsymbol{w}\end{bmatrix}$ is also a multivariate Gaussian with mean
$\boldsymbol{0}$ and inverse-covariance
$\begin{bmatrix}\mathbb{Q}_{\boldsymbol{\phi}}^{-1}&\boldsymbol{0}\\\
\boldsymbol{0}&\sigma^{-2}\mathbb{I}\end{bmatrix}$. Let
$\mathbb{Q}_{\boldsymbol{\phi}}^{-1}=\sigma^{-2}\mathbb{D}$ for some matrix
$\mathbb{D}$. The change-of-variables formula for probability distributions
then tells us that

$\displaystyle\mathbb{\Gamma}_{\boldsymbol{\theta}}=\begin{bmatrix}\mathbb{I}&\mathbb{0}\\\
\mathbb{H}&\mathbb{L}\end{bmatrix}^{\top}\begin{bmatrix}\mathbb{Q}_{\boldsymbol{\phi}}^{-1}&\boldsymbol{0}\\\
\boldsymbol{0}&\kappa^{-1}\mathbb{I}\end{bmatrix}\begin{bmatrix}\mathbb{I}&\mathbb{0}\\\
\mathbb{H}&\mathbb{L}\end{bmatrix}=\frac{1}{\kappa}\begin{bmatrix}\mathbb{I}&\mathbb{0}\\\
\mathbb{H}&\mathbb{L}\end{bmatrix}^{\top}\begin{bmatrix}\mathbb{D}&\boldsymbol{0}\\\
\boldsymbol{0}&\mathbb{I}\end{bmatrix}\begin{bmatrix}\mathbb{I}&\mathbb{0}\\\
\mathbb{H}&\mathbb{L}\end{bmatrix}.$ (118)

(Note that change-of-variables also tells us that
$\log\det\boldsymbol{\Gamma}_{\boldsymbol{\theta}}=\log\det\mathbb{D}-T\log\kappa$,
which we need to compute the term $\boldsymbol{c}_{\theta}$.)

galbraith1974inverses show that for a _stationary_ AR process, we have
$\mathbb{D}=\mathbb{L}_{P}^{\top}\mathbb{L}_{P}-\mathbb{H}_{P}^{\top}\mathbb{H}_{P}$,
where $\mathbb{H}_{P}\in\mathbb{R}^{P\times P}$ is the first $P$ rows of
$\mathbb{H}$ and $\mathbb{L}_{P}\in\mathbb{R}^{P\times P}$ is the top-left
$P\times P$ block of $\mathbb{L}$.

Next, let $\mathbb{D}$ be factorized as
$\mathbb{R}\mathbb{R}^{\top}=\mathbb{D}$, which we can obtain through Cholesky
decomposition. This implies that
$\mathbb{\Gamma}_{\boldsymbol{\theta}}=\mathbb{X}^{\top}\mathbb{X}$, where

$\displaystyle\mathbb{X}:=\frac{1}{\sqrt{\kappa}}\begin{bmatrix}\mathbb{R}&\mathbb{0}\\\
\mathbb{H}&\mathbb{L}\end{bmatrix}.$ (119)

We can verify that $\mathbb{X}$ is a lower triangular and banded matrix with
$P+1$ non-zero bands below the diagonal. Similarly $\mathbb{X}^{\top}$ is an
upper triangular and banded matrix with $P+1$ non-zero bands above the
diagonal. It follows that $\mathbb{\Gamma}_{\boldsymbol{\theta}}$ has $2P+1$
non-zero bands. In turn, this implies that $\mathbb{A}_{\boldsymbol{\theta}}$
in (6) also has $2P+1$ non-zero bands (since the other part of the sum is a
diagonal matrix).

####  F.1.2 Kalman Smoother Implementation of Exact-Gradient EM

Observe that one can easily write (26) as a state-space model by defining the
state

$\displaystyle\boldsymbol{s}_{d}:=\begin{bmatrix}z_{d}\\\ \vdots\\\
z_{d+P-1}\end{bmatrix}$ (120)

for $d=1,\ldots,D$. Then, the noisy AR model of (26) is equivalent to the
following state-space model:

$\displaystyle\boldsymbol{s}_{1}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{0},\mathbb{Q}_{\boldsymbol{\phi}})$
(121) $\displaystyle\boldsymbol{s}_{d}$
$\displaystyle=\mathbb{F}\boldsymbol{s}_{d-1}+\boldsymbol{v}_{d},\quad\boldsymbol{v}_{d}\sim\mathcal{N}(\boldsymbol{0},\mathbb{V})$
(122) $\displaystyle y_{d}$
$\displaystyle=\boldsymbol{c}^{\top}\boldsymbol{s}_{d}+w_{d},\quad
w_{d}\sim\mathcal{N}(0,W)$ (123)

for

$\displaystyle\mathbb{F}:=\begin{bmatrix}0&1&0&\ldots&0\\\ 0&0&1&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\ 0&0&0&\ldots&1\\\
\phi_{P}&\phi_{P-1}&\phi_{P-2}&\ldots&\phi_{1}\end{bmatrix},\quad\mathbb{V}:=\begin{bmatrix}0&0&0&\ldots&0\\\
0&0&0&\ldots&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\ 0&0&0&\ldots&0\\\
0&0&0&\ldots&\kappa\\\ \end{bmatrix},\quad\boldsymbol{c}:=\begin{bmatrix}1\\\
0\\\ 0\\\ \vdots\\\ 0\end{bmatrix},\quad W=\lambda.$ (124)

The Kalman smoother can be used to obtain the distributions
$p(\boldsymbol{s}_{d}|\boldsymbol{y})$ for all $d$, which we can then convert
into our posterior of interest $p(z_{d}|\boldsymbol{y})$. Kalman will require
$D$ steps with complexity $\mathcal{O}(P^{3})$ each because the state size is
$P$. We then follow Bishop & Nasrabadi (2006), Chapter 13.3.2 to compute the
EM objective (5) using the Kalman smoother outputs.

####  F.1.3 Experimental Settings

It is not straightforward to directly optimize $\boldsymbol{\phi}$ over the
space of _stationary_ noisy AR processes. Thus, we parameterize
$\boldsymbol{\phi}:=f(\boldsymbol{\gamma})$, where
$\boldsymbol{\gamma}\in[-1,1]^{P}$ are partial auto-correlations and $f$ is
the transformation defined by barndorff1973parametrization. Gradient descent
is performed over $\boldsymbol{\gamma},\log\kappa$ and $\log\lambda$. For each
algorithm (i.e. gradient EM, probabilistic unrolling), we perform 200
iterations of gradient descent with the Adam optimizer and learning rate
$0.1$.

Each ground-truth value $\gamma^{\star}_{p}$ is randomly initialized between
$[-1,1]$. Similarly, $\kappa^{\star}$ and $\lambda^{\star}$ are randomly
initialized between $[0.1,10]$ (in log-space).

####  F.1.4 Comparison with VAEs

To provide another baseline of comparison for probabilistic unrolling, we fit
the model in (26) using a variational auto-encoder (VAE) (Kingma & Welling,
2013). The decoder of the VAE is the generative model in (26), and the encoder
of the VAE is a black-box deep neural network. To tune the VAE, we search over
different architectures (i.e. 1, 2, 3, 4 layer models), different activations
(i.e. ReLU, Sigmoid, Tanh, LeakyReLU) and adjust weights for the different
parts of the VAE loss (i.e. the weight $\beta$ in $\beta$-VAE
(higgins2017beta)).

However, we find that VAEs perform poorly in parameter recovery for the noisy
AR model. In Table 2, we observe that probabilistic unrolling (like EM)
estimates all the true parameters
$\\{\phi^{\star},\kappa^{\star},\lambda^{\star}\\}$ to within 1% error for
$D=30{,}000$ time points. On the other hand, the (optimally-tuned) VAE obtains
$41.5\pm 14.1\%$ error for $\phi$, $215.93\pm 223.5\%$ error for $\kappa$, and
$37.9\pm 23.75\%$ error for $\lambda$. We hypothesize that this poor
performance is due to the VAE’s biased objective in comparison to EM (and PU),
which suffers from mean-field’s inability to model covariance in the latent
posterior; this is especially detrimental for this problem, because there is a
lot of rich covariance structure across time.

###  F.2 Bayesian Compressed Sensing of Sparse Signals

####  F.2.1 Woodbury Matrix Identity for Exact-Gradient EM

The Woodbury matrix identity is a property from linear algebra that allows us
to compute $\mathbb{\Sigma}_{\boldsymbol{\theta}}$ for Bayesian compressed
sensing by inverting a $M\times M$ matrix instead of a $D\times D$ one. Since
typically $M<D$ in compressed sensing, this can lead to computational savings
in time from a practical standpoint (but perhaps not an asymptotic one,
because $M$ often needs to grow linearly with $D$ in compressed sensing
applications). Note that it does not lead to computational savings in space,
because the full matrix $\mathbb{\Sigma}_{\boldsymbol{\theta}}$ is still
computed in the end.

Instead of computing $\mathbb{\Sigma}_{\boldsymbol{\theta}}$ through (7), we
equivalently have

$\displaystyle\mathbb{\Sigma}_{\boldsymbol{\theta}}=\mathbb{\Gamma}_{\boldsymbol{\theta}}^{-1}-\mathbb{\Gamma}_{\boldsymbol{\theta}}^{-1}\mathbb{\Phi}_{\boldsymbol{\theta}}^{\top}\mathbb{\Omega}^{\top}(\mathbb{\Omega}\mathbb{\Psi}_{\boldsymbol{\theta}}^{-1}\mathbb{\Omega}^{\top}+\mathbb{\Omega}\mathbb{\Phi}_{\boldsymbol{\theta}}\mathbb{\Gamma}_{\boldsymbol{\theta}}^{-1}\mathbb{\Phi}_{\boldsymbol{\theta}}^{\top}\mathbb{\Omega}^{\top})^{-1}\mathbb{\Omega}\mathbb{\Phi}_{\boldsymbol{\theta}}\mathbb{\Gamma}_{\boldsymbol{\theta}}^{-1}.$
(125)

Note that for Bayesian compressed sensing, both
$\mathbb{\Gamma}_{\boldsymbol{\theta}}^{-1}$ and
$\mathbb{\Psi}_{\boldsymbol{\theta}}^{-1}$ can be cheaply computed, because
these are both diagonal matrices.

####  F.2.2 Experimental Settings

The NIST dataset is accessed at https://www.nist.gov/srd/nist-special-
database-19.

We scale all raw image pixels in NIST from $[0,255]$ to the range $[0,1]$. We
add independent, pixel-wise Gaussian noise to the undersampled 2D Fourier
transform with standard deviation $\sigma=0.01$. During model fitting, we use
the Adam optimizer with learning rate 1.0. We optimize the parameters
$\boldsymbol{\alpha},\beta$ in log-space. Each component of
$\log\boldsymbol{\alpha}$ is initialized as randomly drawn from
$\mathcal{N}(0,1)$. The value $\log\beta$ is initialized as 0.

For probabilistic unrolling, we use the preconditioned conjugate gradient
algorithm. The preconditioner we employ is the diagonal preconditioner
$\mathbb{M}$ introduced in Lin et al. (2022b) for sparse Bayesian learning;
$\mathbb{M}^{-1}$ is a diagonal matrix with diagonal $\boldsymbol{m}$, where

$\displaystyle m_{j}:=\frac{1}{\alpha_{j}+\beta}.$ (126)

####  F.2.3 Sample Images

In Figure 3, we provide sample images of the true signal
$\boldsymbol{\tilde{z}}^{(n)}$, its 15% undersampled and noisy Fourier
transform measurement $\boldsymbol{\tilde{y}}^{(n)}$, and a reconstruction
provided by probabilistic unrolling $\boldsymbol{\mu}^{(n)}$.

Figure 3: From left to right: the true signal, the Fourier measurement, and
the probabilistic unrolling reconstruction.

####  F.2.4 Results Breakdown by Digit Type

Table 5 presents a breakdown of the compressed sensing results by digit type.
The average is given in Table 3.

Table 5: Compressed sensing results broken down by digit type. Digit Type | $r(\boldsymbol{\mu}^{\text{EM}},\boldsymbol{\tilde{z}})$ | $r(\boldsymbol{\mu}^{\text{PU}},\boldsymbol{\tilde{z}})$ | EM Time | PU Time
---|---|---|---|---
0 | 3.89% | 4.01 % | 1475 s | 21 s
1 | 4.76% | 4.31 % | 1483 s | 21 s
2 | 4.61% | 4.43 % | 1473 s | 21 s
3 | 3.96% | 3.90 % | 1508 s | 21 s
4 | 7.20% | 8.56 % | 1485 s | 21 s
5 | 5.22% | 4.88 % | 1494 s | 21 s
6 | 4.17% | 3.94 % | 1510 s | 21 s
7 | 4.67% | 4.59 % | 1467 s | 21 s
8 | 4.27% | 4.08 % | 1468 s | 21 s
9 | 4.83% | 4.46 % | 1448 s | 21 s

Figure 4: Gradient convergence as a function of iterations $I$.

####  F.2.5 Super-Efficiency of Network Gradient

In Figure 4, we empirically show that the network gradient (17) converges
faster to the Monte Carlo gradient (10) than the output gradient (16) for our
preconditioned conjugate gradient solver.

###  F.3 Collaborative Filtering through Factor Analysis

####  F.3.1 Experimental Settings

The MovieLens datasets are accessed at
https://grouplens.org/datasets/movielens/.

We follow the general experimental setup of Sedhain et al. (2015). For each
dataset, we perform a 90%-10% train-test split of the ratings data. We then
set aside 10% of the training set as a validation set. Both models (gradient
EM and probabilistic unrolling) are trained with the Adam optimizer, learning
rate 0.001, and a user mini-batch size of 25 (with gradient accumulation over
four mini-batches for an overall batch size of 100). For the ML-1M dataset, we
train the model for 20 epochs and evaluate the model on the validation set
every 500 gradient steps. For the ML-10 M and ML-25 M datasets, we train the
model for 5 epochs and evaluate the model on the validation set every 2,000
gradient steps. The checkpoint with the best validation set RMSE is used for
final evaluation on the test set.

During training, we only train on users with at least one rating in the
training set. Thus, there may be users in the validation/test sets that do not
appear during training. Following Sedhain et al. (2015), we always predict a
rating of 3 for these users. All other model rating predictions are clipped to
be in the range [1.0, 5.0] before evaluation of RMSE.

####  F.3.2 Comparison with VAEs

Similar to Appendix F.1.4, we report results for a VAE baseline on the
collaborative filtering task. Searching over architectural options, we find
that the best-performing architecture had a two-layer encoder with $3{,}000$
units in the hidden layer and ReLU activations. However, even with our
extensive tuning, we observe that VAEs generally perform worse and have higher
computational costs in comparison to probabilistic unrolling. The VAE obtains
0.8849 RMSE (compared to 0.8436 for PU) on MovieLens 1-m and 0.8366 RMSE
(compared to 0.7796 for PU) on MovieLens 10-m. The VAE also has approximately
1.5x the time cost and 2x the memory cost of PU (due to the separate inference
encoder network).


The Extended Abstracts of
The 4th Seminar on Functional Analysis and its Applications
2-3rd March 2016, Ferdowsi University of Mashhad, Iran

# Relative Reproducing Kernels in Vector-Valued Hilbert and Banach Spaces

Ali Ebadian1,∗00footnotetext: ∗ Speaker.Saeed Hashemi Sababe2 1 Department of
Mathematics, Payame Noor University (PNU), Iran;
<EMAIL_ADDRESS>
2 Department of Mathematics, Payame Noor University (PNU), Iran;
Hashemi<EMAIL_ADDRESS>
###### Abstract.

This paper is devoted to the study of vector valued reproducing kernel Hilbert
spaces. We focus on reproducing kernels in vector-valued reproducing kernel
Hilbert spaces. In particular we extend reproducing kernels to relative
reproducing kernels and prove some theorems in this subject.

###### Key words and phrases:

Hilbert space, reproducing kernels, semi-normed spaces.

###### 2010 Mathematics Subject Classification:

Primary 47B32; Secondary 47A70.

##  1\. Introduction

The purpose of this paper is to establish the notion of vector-valued relative
reproducing kernels in Hilbert and Banach spaces. In learning theory,
reproducing kernel Hilbert spaces (RKHS) and reproducing kernel Banach spaces
(RKBS) are important tools for designing learning algorithms. The basic object
is a Hilbert or Banach space of functions $f$ from a set $X$ into a normed
vector space $\mathcal{Y}$ with the property that, for any $x\in X,|f(x)|\leq
C_{x}|f|$ for a positive constant $C_{x}$ independent of $f$. In continues, we
briefly review basic definitions. More details can be found in [1, 2, 3, 4,
5].

Given a set $X$ and a normed vector space $\mathcal{Y}$, a map $K:X\times
X\rightarrow\mathcal{L}(\mathcal{Y})$ is called a $\mathcal{Y}$-reproducing
kernel if

$\sum_{i,j=1}^{n}\langle K(x_{i},x_{j})y_{j},y_{i}\rangle\geq 0$ (1.1)

or any $x_{1},\dots,x_{n}$ in $X$, $y_{1},\dots,y_{n}$ in $\mathcal{Y}$ and
$n\geq 1$. Given $x\in X$,
$K_{x}:\mathcal{Y}\rightarrow\mathcal{F}(X,\mathcal{Y})$ denotes the linear
operator whose action on a vector $y\in\mathcal{Y}$ is the function
$K_{xy}\in\mathcal{F}(X,\mathcal{Y})$ defined by

$(K_{xy})(t)=K(t,x)y\qquad t\in X.$ (1.2)

Given a $\mathcal{Y}$-reproducing kernel $K$, there is a unique Hilbert space
$\mathcal{H}_{K}\subset\mathcal{F}(X,\mathcal{Y})$ satisfying

$\displaystyle K_{x}\in L(\mathcal{Y},\mathcal{H}_{K})\qquad x\in X$ (1.3)
$\displaystyle f(x)=K_{x}^{∗}f\qquad x\in X,f\in\mathcal{H}_{K},$ (1.4)

where $K_{x}^{*}:\mathcal{H}_{K}\rightarrow\mathcal{Y}$ is the adjoint of
$K_{x}$. The space $\mathcal{H}_{K}$ is called the reproducing kernel Hilbert
space associated with $K$, the corresponding scalar product and norm are
denoted by $\langle.,.\rangle_{K}$ and $\|.\|_{K}$, respectively. As a
consequence of $(\ref{3})$, we have that

$\displaystyle K(x,t)=K_{x}^{∗}K_{t}\qquad x,t\in X$ (1.5)
$\displaystyle\mathcal{H}_{K}=\overline{span}\\{K_{xy}|x\in X,y\in Y\\}.$
(1.6)

Let $\mathcal{H}$ be a Hilbert space of functions defined on the set $X$. We
say that it is a relative reproducing kernel Hilbert space if there exists a
function $M_{x,y}$ from $X\times X$ into $\mathcal{H}$ such that

$F(x)-F(y)=\langle F,M_{x,y},\mathcal{H}\rangle\qquad\forall x,y\in X,\forall
F\in\mathcal{H}.$ (1.7)

A semi-inner product on a Banach space $V$ is a function from $V\times V$ to
$\mathbb{C}$, denoted by $[.,.]_{V}$, such that for all $u,v,w\in V$ and
$\alpha,\beta\in\mathbb{C}$

  1. (1)

(linearity ) $[\alpha f+\beta g,h]_{V}=\alpha[f,h]_{V}+\beta[g,h]_{V}$;

  2. (2)

(positivity) $[f,f]_{V}>0$ for $f\neq 0$;

  3. (3)

(conjugate homogeneity) $[f,\alpha g]_{V}=\alpha[f,g]_{V}$;

  4. (4)

(Cauchy-Schwartz inequality) $|[f,g]_{V}|\leq[f,f]_{V}^{1/2}[g,g]_{V}^{1/2}$.

A semi-inner product $[.,.]_{V}$ on $V$ is said to be compatible if

$[f,f]_{V}^{1/2}=\|f\|_{V}\qquad\forall f\in V,$ (1.8)

where $\|.\|_{V}$ denotes the norm on $V$. Every Banach space has a compatible
semi-inner product. Let $[.,.]_{V}$ be a compatible semi-inner product on $V$.
Then one sees by the Cauchy-Schwartz inequality that for each
$f\in\mathcal{B}$, the linear functional $f^{*}$ on $V$ defined by

$f^{*}(g):=[g,f]_{V},\qquad g\in V$ (1.9)

is bounded on $V$. In other words, $f^{*}$ lies in the dual space
$\mathcal{B}^{∗}$ of $\mathcal{B}$. Moreover, we have

$\|f^{*}\|_{V^{*}}=\|f\|_{V}$ (1.10)

Introduce the duality mapping $\mathcal{I}_{V}$ from $V$ to $V^{∗}$ by setting

$\mathcal{I}_{V}(f):=f^{∗},\qquad f\in V.$ (1.11)

Let $\Lambda$ be a Banach space. A space $\mathcal{B}$ is called a Banach
space of $\Lambda$-valued functions on $X$ if it consists of certain functions
from $X$ to $\Lambda$ and the norm on $\mathcal{B}$ is compatible with point
evaluations in the sense that

$\|f\|_{B}=0\Longleftrightarrow f(x)=0\quad\forall x\in X.$ (1.12)

We call $\mathcal{B}$ a $\Lambda$-valued RKBS on $X$ if both $\mathcal{B}$ and
$\Lambda$ are uniform and $\mathcal{B}$ is a Banach space of functions from
$X$ to $\Lambda$ such that for every $x\in X$, the point evaluation
$\delta_{x}:\mathcal{B}\rightarrow\lambda$ defined by

$\delta_{x}(f):=f(x),\qquad f\in\mathcal{B},$ (1.13)

is continuous from $\mathcal{B}$ to $\Lambda$.

##  2\. Main Results

######  Definition 2.1.

Let $X$ be set an arbitrary set, $\mathcal{Y}$ a normed vector space and
$\mathcal{H}$ a vector-valued reproducing kernel Hilbert space with $K$ as its
reproducing kernel satisfying in properties $(\ref{1},\ref{2})$. We say that
it is a vector-valued relative reproducing kernel Hilbert space (RRKHS) if
there exists a function $M_{x,y}$ such that

$\displaystyle M_{x,y}\in\mathcal{L}(\mathcal{Y},\mathcal{H})\qquad\forall
x,y\in X$ (2.1) $\displaystyle f(y)-f(x)=M_{x,y}^{*}f\qquad\forall x,y\in
X,\forall f\in\mathcal{H}$ (2.2)

where $M_{x,y}^{*}:\mathcal{H}\rightarrow\mathcal{Y}$ is the adjoint of
$M_{x,y}$.

######  Definition 2.2.

Let $\mathcal{B}$ be a Banach space of $\Lambda$-valued functions on $X$. We
call $\mathcal{B}$ a $\Lambda$-valued relative reproducing kernel Banach space
(RRKBS) on $X$ if both $\mathcal{B}$ and $\Lambda$ are uniform and
$\mathcal{B}$ is a Banach space of functions from $X$ to $\Lambda$ such that
for every $x\in X$, the function $\zeta_{x,y}:\mathcal{B}\rightarrow\lambda$
defined by

$\zeta_{x,y}(f):=f(y)-f(x),\qquad\forall f\in\mathcal{B},$ (2.3)

is continuous from $\mathcal{B}$ to $\Lambda$.

######  Theorem 2.3.

The function $M_{x,y}$ in (2.1) is unique and satisfies

$M_{x_{1},x_{2}}(t)+M_{x_{2},x_{3}}(t)=M_{x_{1},x_{3}}(t)\qquad\forall
x_{1},x_{2},x_{3},t\in X.$ (2.4)

###### Proof.

The uniqueness of the function $M_{x_{1},x_{2}}$ follows from Riesz’
representation theorem. This uniqueness and following equality imply (2.4).

$f(x_{1})-f(x_{3})=f(x_{1})-f(x_{2})+f(x_{2})-f(x_{3})$ (2.5)

∎

######  Theorem 2.4.

Previews theorem also holds in RRKBS but the function is not necessarily
unique in these spaces.

###### Proof.

similarly to theorem 2.3. ∎

Obvously the Riesz’ representation theorem does not hold in general Banach
spaces but with adding some conditions, this theorem can be hold. So in
general case, the above uniqueness does not hold.

######  Theorem 2.5.

Let $\mathcal{B}$ be a $RRKBS$ on a locally compact Hausdorf Banach space
$\Lambda$. Then the function defined in (2.3) is unique.

###### Proof.

In this case, a version of Riesz’ representation theorem holds and it implies
the uniqueness. see Theorem 1 in [3]. ∎

######  Theorem 2.6.

Let $\mathcal{H}$ be a Hilbert space and there exist functions $K$ and $M$
with reproducing and relative reproducing property, respectively, then

$\mathcal{H}_{M}\subseteq\mathcal{H}_{K}\subseteq\mathcal{H}$ (2.6)

## Acknowledgement

Research supported in part by Kavosh Alborz Institute of Mathematics and
Applied Sciences and was accomplished while the second author was visiting the
University of Iowa.

## References

  * [1] D. Alpay , P. Jorgensen A N D D. Volok, Relative reproducing kernel Hilbert spaces, Proceedings of the American mathematical society Volume 142, Number 11, November 2014, Pages 3889 – 3895.
  * [2] A. Ebadian, S. Hashemi Sababe and Sh. Najafzadeh, On generalized reproducing kernel Hilbert spaces , Submited in Operators and Matrices, arXiv:1504.04668.
  * [3] A. Ebadian, S. Hashemi Sababe and Sh. Najafzadeh, Generalized reproducing kernel Banach spaces in machine learning, Submited in Bulletin of the Korean Mathematical Society.
  * [4] S. Saitoh, Theory of reproducing kernels and its applications,volume 189 of Pitman Research Notes in Mathematics Series. Longman Scientific & Technical, Harlow, 1988.
  * [5] Haizhang Zhang and Jun Zhang, Vector-valued Reproducing Kernel Banach Spaces with Applications to Multi-task Learning, arXiv:1111.1037v2.