# Linear Canonical Transform On Boehmian Space

S. K. Panchal and Pravinkumar V. Dole
Department of Mathematics,
Dr. Babasaheb Ambedkar Marathwada University,
Aurangabad-431004 (M.S.) India.
E-mail ID-<EMAIL_ADDRESS><EMAIL_ADDRESS>
Abstract: The aim of this paper is to constructs Boehmian space, the linear
canonical transform for Boehmians is define and to study its properties.

AMS Subject Classification: 44A35, 44A40, 46F12, 46F99.

Key Words: Linear canonical transform, Convolution, Distributions, Boehmians.

##  1 Introduction

The most recent generalizations of functions is the theory of Boehmians. The
idea of construction of Boehmians was initiated by the concept of regular
operators introduced by Boehme [1]. Regular operators form a subalgebra of the
field of Mikusinski operators and they include only such functions whose
support is bounded from the left. In a concrete case, the space of Boehmians
contains all regular operators, all distributions and some objects which are
neither operators nor distributions. The space of Boehmians is the new class
of generalized functions which opened the new door to area of research in
mathematics. The construction of Boehmians is given by Mikusinski and
Mikusinski [5, 6, 8, 9]. Mikusinski and Nemzer has studied Fourier and Laplace
transform for Boehmians in [7] and [11] respectively. Zayed [15] extended the
fractional Fourier transform to class of integrable Boehmians. Singh studied
fractional integrals of fractional Fourier transform for integrable Boehmians
in [13]. The Fourier, Laplace and fractional Fourier transforms are the
special cases of linear canonical transform (LCT) and has many applications in
several areas, like signal processing and optics [2]. This lead to study
linear canonical transform for integrable Boehmians in [4]. In one hand we
constructs Boehmian spaces and other hand defined linear canonical transform
for Boehmians. Further, we obtain its properties like one-to-one, onto,
continuous from one Boehmian space to another Boehmian space and other basic
properties in the space of Boehmians.

The linear canonical transform of real valued function $f$ is defined [3, 10]
as,

$\displaystyle\mathcal{L}_{A}[f(t)](u)=F_{A}(u)=\left\\{{\begin{array}[]{*{20}{l}}{\sqrt{\frac{1}{2\pi
ib}}\int_{-\infty}^{\infty}e^{\frac{i}{2}[\frac{a}{b}t^{2}-\frac{2}{b}ut+(\frac{d}{b})u^{2}]}f(t)dt\,\,for\,\,b\neq
0,}\\\ {\sqrt{d}e^{\frac{i}{2}cdu^{2}}f(du)\qquad\qquad\qquad\qquad
for\,\,b=0,}\end{array}}\right.$ (1.3)

where $\mathcal{L}_{A}$ is the unitary linear canonical transform operator
with parameter $A=(a,b,c,d)$, $a,b,c,d$ are real number satisfying $ad-bc=1$.
The inverse transform for linear canonical transform is given by a linear
canonical transform having the parameter $A^{-1}=(d,-b,-c,a)$ and
$\mathcal{L}_{A^{-1}}$ is the inverse LCT operator. For value of parameters
as, $a=cos\theta,b=sin\theta,c=-sin\theta,d=cos\theta$ then LCT become
fractional Fourier transform, in particular, when $\theta=\frac{\pi}{2}$ then
LCT become Fourier transform and $a=0,b=i,c=i,d=0$ then LCT becomes Laplace
transform.

Let $\mathcal{L}^{1}(\mathbb{R})$ be the space of all complex valued
absolutely integrable functions on $\mathbb{R}$ with norm
$||f||_{1}=\int_{\mathbb{R}}|f(t)|dt\leq M_{1}$ and
$\mathcal{L}^{2}(\mathbb{R})$ be the space of all complex valued absolutely
square integrable functions on $\mathbb{R}$ with norm
$||g||_{2}=\big{(}\int_{\mathbb{R}}|g(t)|^{2}dt\big{)}^{\frac{1}{2}}\leq
M_{2}$, for some $M_{1},M_{2}>0$. Let
$\mathcal{L}^{1}(\mathbb{R})\cap\mathcal{L}^{2}(\mathbb{R})$ is denoted by
$\mathcal{L}^{1,2}(\mathbb{R})$.

###### Definition 1.1

[12](Regular Distributions) Let $f$ be the locally integrable function, i.e.
absolutely integrable on every finite interval on $\mathbb{R}$, then
distribution generated by $f$ is called regular distributions.

We see that $f\in\mathcal{L}^{1,2}(\mathbb{R})$ then $\mathcal{L}_{A}(f)$ and
$\mathcal{L}_{A^{-1}}(f)$ are the members of $\mathcal{L}^{1,2}(\mathbb{R})$.

###### Definition 1.2

[3] Let the weight function $W(t,\tau)=e^{i\tau(\tau-t)\frac{a}{b}}$. For any
two function $f$ and $g$ the convolution operation $*^{A}$ is defined as,

$\displaystyle
h(t)=(f*^{A}g)(t)=\int_{-\infty}^{\infty}f(\tau)g(t-\tau)W(t,\tau)d\tau$ (1.4)

###### Theorem 1.1

[3] (New Convolution Theorem)
Let $h(t)=(f*^{A}g)(t)$ and $H_{A}(u),F_{A}(u),G_{A}(u)$ denote the linear
canonical transform of $h(t),f(t)$ and $g(t)$ respectively, then

$\displaystyle H_{A}(u)=\sqrt{2i\pi
b}\,e^{-i(\frac{du^{2}}{2b})}F_{A}(u)G_{A}(u).$ (1.5)

##  2 Preliminary Results

In this section we obtain some results which are require to construct the
Boehmian space.

###### Lemma 2.1

Let $f\in\mathcal{L}^{1}(\mathbb{R})$ and $g\in\mathcal{L}^{2}(\mathbb{R})$
then the $(f*^{A}g)$ is in $\mathcal{L}^{2}(\mathbb{R})$.

###### Lemma 2.2

The space $(\mathcal{L}^{1,2}(\mathbb{R}),*^{A})$ is commutative semi group.

###### Theorem 2.1

(Plancherel type theorem) Let the sequence
$f_{n}\in\mathcal{L}^{1,2}(\mathbb{R})$ and $f_{n}\rightarrow f$ on
$\mathcal{L}^{2}(\mathbb{R})$ then
$\mathcal{L}_{A}(f_{n})\rightarrow\mathcal{L}_{A}(f)$ in
$\mathcal{L}^{2}(\mathbb{R})$ as $n\rightarrow\infty$.

###### Definition 2.1

Analogous to Plancherel type theorem for $f\in\mathcal{L}^{2}(\mathbb{R})$, we
define $\mathcal{L}_{A}(f)$ by
$\mathcal{L}^{2}-\lim_{n\rightarrow\infty}\mathcal{L}_{A}(f_{n})$, where
$f_{n}\in\mathcal{L}^{1,2}(\mathbb{R})$.

Let $\bigtriangledown$ be the set of all sequences of continuous real
functions $\\{\delta_{n}\\}$ from $\mathcal{L}^{1,2}(\mathbb{R})$ having
compact support on $\mathbb{R}$ with the following properties:

  1. (i)

$\quad\int_{\mathbb{R}}e^{i\frac{at^{2}}{2b}}\delta_{n}(t)dt=1$,
$\forall\,n\in\mathbb{N}$,

  2. (ii)

$\quad\lim_{n\rightarrow\infty}\int_{|t|>\epsilon}|\delta_{n}(t)|dt=0$ for
each $\epsilon>0$.

The members of $\bigtriangledown$ are called delta sequences.

###### Example 2.1

Let $a,b\in\mathbb{R};b\neq 0$, consider the sequence

$\displaystyle\delta_{n}(t)=\left\\{{\begin{array}[]{*{20}{l}}{e^{-i\frac{at^{2}}{2b}}t\qquad\qquad\qquad
for\quad 0\leq t\leq\frac{1}{n},}\\\
{e^{-i\frac{at^{2}}{2b}}n^{2}(\frac{2}{n}-t)\qquad\quad\,\,for\quad\frac{1}{n}\leq
t\leq\frac{2}{n},}\\\ {0\qquad\qquad\qquad\qquad\qquad\quad
otherwise}.\end{array}}\right.$

###### Lemma 2.3

Let $\\{\varphi_{n}\\},\\{\psi_{n}\\}\in\bigtriangledown$ then
$(\varphi_{n}*^{A}\psi_{n})\in\bigtriangledown$ for all $n\in\mathbb{N}$.

###### Lemma 2.4

Let $f\in\mathcal{L}^{1,2}(\mathbb{R})$ and
$\\{\psi_{n}\\}\in\bigtriangledown$ then $f*^{A}\psi_{n}\rightarrow f$ as
$n\rightarrow\infty$ in $\mathcal{L}^{2}(\mathbb{R})$.

##  3 LCT For Boehmians

A pair of sequences $(f_{n},\varphi_{n})$ is called a quotient of the
sequences, denoted by $f_{n}/\varphi_{n}$, where each $n\in\mathbb{N}$,
$f_{n}\in\mathcal{L}^{1,2}(\mathbb{R})$ and
$\\{\varphi_{n}\\}\in\bigtriangledown$ such that
$f_{m}*^{A}\varphi_{n}=f_{n}*^{A}\varphi_{m}$ holds
$\forall\,m,n\in\mathbb{N}$. Two quotients of sequences $f_{n}/\varphi_{n}$
and $g_{n}/\psi_{n}$ are equivalent if
$f_{n}*^{A}\psi_{n}=g_{n}*^{A}\varphi_{n}$ $\forall\,n\in\mathbb{N}$. This is
an equivalence relation. The equivalence class of quotient of sequences is
called a Boehmian. The space of all Boehmians is denoted by
$\mathcal{B}_{\mathcal{L}^{1,2}}=\mathcal{B}_{\mathcal{L}^{1,2}}(\mathcal{L}^{1,2}(\mathbb{R}),\bigtriangledown,*^{A})$
and the members of $\mathcal{B}_{\mathcal{L}^{1,2}}$ are denoted by
$F=[f_{n}/\varphi_{n}]$. The function $f\in\mathcal{L}^{1,2}(\mathbb{R})$ can
be identified with the Boehmian $[(f*^{A}\delta_{n})/\delta_{n}]$, where
$\\{\delta_{n}\\}$ is the delta sequence. Let $F=[f_{n}/\varphi_{n}]$, then
$F*^{A}\delta_{n}=f_{n}\in\mathcal{L}^{1,2}(\mathbb{R})$
$\forall\,n\in\mathbb{N}$.

###### Definition 3.1

A sequence of Boehmians $F_{n}$ is called $\Delta-$convergent to a Boehmian
$F$ ($\Delta-\lim F_{n}=F$) if there exist a delta sequence $\\{\delta_{n}\\}$
such that $(F_{n}-F)*^{A}\delta_{n}\in\mathcal{L}^{1,2}(\mathbb{R})$, for
every $n\in\mathbb{N}$ and that $\|(F_{n}-F)*^{A}\delta_{n}\|_{2}\rightarrow
0$ as $n\rightarrow\infty$.

###### Definition 3.2

A sequence of Boehmians $F_{n}$ is called $\delta-$convergent to a Boehmian
$F$ ($\delta-\lim F_{n}=F$) if there exist a delta sequence $\\{\delta_{n}\\}$
such that $F_{n}*^{A}\delta_{k}\in\mathcal{L}^{1,2}(\mathbb{R})$ and
$F*^{A}\delta_{k}\in\mathcal{L}^{1,2}(\mathbb{R})$ for every
$n,k\in\mathbb{N}$ and that $\|(F_{n}-F)*^{A}\delta_{k}\|_{2}\rightarrow 0$ as
$n\rightarrow\infty$ for each $k\in\mathbb{N}$.

Let $\\{\delta_{n}\\}$ is a delta sequence, then $\delta_{n}/\delta_{n}$
represents an Boehmian. Since the Boehmian $[\delta_{n}/\delta_{n}]$
corresponds to Dirac delta distribution $\delta$, all the derivative of
$\delta$ are also Boehmian. If $\\{\delta_{n}\\}$ is infinitely differentiable
and bounded support, then the $k^{th}$ derivative of $\delta$ is define by
$\delta^{(k)}=[\delta_{n}^{(k)}/\delta_{n}]\in\mathcal{B}_{\mathcal{L}^{1,2}}$,
for each $k\in\mathbb{N}$. The $k^{th}$ derivative of Boehmian
$F\in\mathcal{B}_{\mathcal{L}^{1,2}}$ is define by
$F^{(k)}=F*^{A}\delta^{(k)}$.The scalar multiplication, addition and
convolution in $\mathcal{B}_{\mathcal{L}^{1,2}}$ are define as,

$\displaystyle\lambda[f_{n}/\varphi_{n}]$ $\displaystyle=[\lambda
f_{n}/\varphi_{n}]$ $\displaystyle[f_{n}/\varphi_{n}]+[g_{n}/\psi_{n}]$
$\displaystyle=[(f_{n}*^{A}\psi_{n}+g_{n}*^{A}\varphi_{n})/\varphi_{n}*^{A}\psi_{n}]$
$\displaystyle[f_{n}/\varphi_{n}]*^{A}[g_{n}/\psi_{n}]$
$\displaystyle=[(f_{n}*^{A}g_{n})/(\varphi_{n}*^{A}\psi_{n})].$

###### Lemma 3.1

Let $\Delta-\lim F_{n}=F$ in $\mathcal{B}_{\mathcal{L}^{1,2}}$ then
$\Delta-\lim F_{n}^{(k)}=F^{(k)}$ for $\forall\,k\in\mathbb{N}$ in
$\mathcal{B}_{\mathcal{L}^{1,2}}$.

Let
$\bigtriangledown_{0}=\\{\mathcal{L}_{A}(\delta_{n});\\{\delta_{n}\\}\in\bigtriangledown\\}$
be the space of complex valued functions on $\mathbb{R}$, the operation
$\cdot$ is pointwise multiplication and $C_{0}(\mathbb{R})$ be the space of
all continuous functions vanishing at infinity on $\mathbb{R}$ then we
construct the another space of Boehmians, denoted by
$\mathcal{B}_{\bigtriangledown}=\mathcal{B}_{\bigtriangledown}(\mathcal{L}^{2}(\mathbb{R}),C_{0}(\mathbb{R})\cap\mathcal{L}^{2}(\mathbb{R}),\cdot,\bigtriangledown_{0})$.
This is the range of linear canonical transform on
$\mathcal{B}_{\mathcal{L}^{1,2}}$ and each element of
$\mathcal{B}_{\bigtriangledown}$ is denoted by
$\mathcal{L}_{A}(f_{n})/\mathcal{L}_{A}(\delta_{n})$ for all $n\in\mathbb{N}$,
where $\\{f_{n}\\}\in\mathcal{L}^{1,2}(\mathbb{R})$.

###### Lemma 3.2

Let $f,g\in\mathcal{L}^{2}(\mathbb{R});\varphi,\psi\in C_{0}(\mathbb{R})$ and
$\lambda\in\mathbb{C}$ then
(i) $f\cdot\varphi\in\mathcal{L}^{2}(\mathbb{R})$
(ii) $(f+g)\cdot\varphi=f\cdot\varphi+f\cdot\varphi$
(iii) $(\lambda f)\cdot\varphi=\alpha(f\cdot\varphi)$
(iv) $f\cdot(\varphi\cdot\psi)=(f\cdot\varphi)\cdot\psi$.

###### Lemma 3.3

Let $f_{n}\rightarrow f$ as $n\rightarrow\infty$ in
$\mathcal{L}^{2}(\mathbb{R})$ and $\varphi\in C_{0}(\mathbb{R})$ then
$f_{n}\cdot\varphi\rightarrow f\cdot\varphi$ in $\mathcal{L}^{2}(\mathbb{R})$.

###### Lemma 3.4

Let $\\{\delta_{n}\\}\in\bigtriangledown$ then $\mathcal{L}_{A}(\delta_{n})$
converges uniformly on each compact set to a constant function $1$ in
$\mathcal{L}^{2}(\mathbb{R})$.

###### Lemma 3.5

Let $f_{n}\longrightarrow f$ as $n\longrightarrow\infty$ in
$\mathcal{L}^{1,2}(\mathbb{R})$ and
$\mathcal{L}_{A}(\varphi_{n})\in\bigtriangledown_{0}$ then
$f_{n}\cdot\mathcal{L}_{A}(\varphi_{n})\rightarrow f$ in
$\mathcal{L}^{2}(\mathbb{R})$.

###### Lemma 3.6

Let
$\mathcal{L}_{A}(\varphi_{n}),\mathcal{L}_{A}(\psi_{n})\in\bigtriangledown_{0}$
then
$\mathcal{L}_{A}(\varphi_{n})\cdot\mathcal{L}_{A}(\psi_{n})\in\bigtriangledown_{0}$.

Proof: Let $\mathcal{L}_{A}(\varphi_{n}),\mathcal{L}_{A}(\psi_{n})\in
C_{0}(\mathbb{R})$ From theorem (1.1) and lemma (2.3) we get
$\mathcal{L}_{A}(\varphi_{n})\cdot\mathcal{L}_{A}(\psi_{n})=\frac{e^{\frac{i}{2}(\frac{d}{b})u^{2}}}{\sqrt{2\pi
ib}}\mathcal{L}_{A}(\varphi_{n}*^{A}\psi_{n})\in\bigtriangledown_{0}$.$\hfill\blacksquare$

###### Definition 3.3

Let $\\{f_{n}\\}\in\mathcal{L}^{1,2}(\mathbb{R})$ and
$\\{\delta_{n}\\}\in\bigtriangledown$, we define the linear canonical
transform
$\mathcal{L}_{A}:\mathcal{B}_{\mathcal{L}^{1,2}}\longrightarrow\mathcal{B}_{\bigtriangledown}$
as

$\displaystyle\mathcal{L}_{A}[f_{n}/\delta_{n}]=\mathcal{L}_{A}(f_{n})/\mathcal{L}_{A}(\delta_{n})\qquad
for\quad[f_{n}/\delta_{n}]\in\mathcal{B}_{\mathcal{L}^{1,2}}.$ (3.1)

The linear canonical transform on $\mathcal{B}_{\mathcal{L}^{1,2}}$ is well
defined. Indeed if $[f_{n}/\delta_{n}]\in\mathcal{B}_{\mathcal{L}^{1,2}}$,
then $f_{n}*^{A}\delta_{m}=f_{m}*^{A}\delta_{n}$ for all $m,n\in\mathbb{N}$.
Applying the linear canonical transform on both sides, we get
$\mathcal{L}_{A}(f_{n})\mathcal{L}_{A}(\delta_{m})=\mathcal{L}_{A}(f_{m})\mathcal{L}_{A}(\delta_{n})$
for all $m,n\in\mathbb{N}$ and hence
$\mathcal{L}_{A}(f_{n})/\mathcal{L}_{A}(\delta_{n})\in\mathcal{B}_{\bigtriangledown}$.
Further if
$[f_{n}/\psi_{n}]=[g_{n}/\delta_{n}]\in\mathcal{B}_{\mathcal{L}^{1,2}}$ then
we have $f_{n}*^{A}\delta_{n}=g_{n}*^{A}\psi_{n}$ for all $n\in\mathbb{N}$.
Again applying the linear canonical transform on both sides, we get
$\mathcal{L}_{A}(f_{n})\mathcal{L}_{A}(\delta_{n})=\mathcal{L}_{A}(g_{n})\mathcal{L}_{A}(\psi_{n})$
for all $n\in\mathbb{N}$. i.e.
$\mathcal{L}_{A}(f_{n})/\mathcal{L}_{A}(\psi_{n})=\mathcal{L}_{A}(g_{n})/\mathcal{L}_{A}(\delta_{n})$
in $\mathcal{B}_{\bigtriangledown}$.

###### Lemma 3.7

Let $[f_{n}/\varphi_{n}]\in B_{\mathcal{L}^{1,2}}$ then the linear canonical
transform of the sequence

$\displaystyle\mathcal{L}_{A}[f_{n}](u)=\sqrt{\frac{1}{2\pi
ib}}e^{\frac{i}{2}(\frac{d}{b})u^{2}}\int_{-\infty}^{\infty}e^{\frac{-i}{b}ut}e^{\frac{i}{2}\frac{a}{b}t^{2}}f_{n}(t)dt$
(3.2)

converges uniformly on each compact set in $\mathbb{R}$.

###### Definition 3.4

In view of the above proof of lemma (3.7), the linear canonical transform of
Boehmian in the space of continuous functions on $\mathbb{R}$ is define as,

$\displaystyle\mathcal{L}_{A}[F]=\lim_{n\rightarrow\infty}\mathcal{L}_{A}(f_{n}).$

###### Theorem 3.1

The linear canonical transform
$\mathcal{L}_{A}:\mathcal{B}_{\mathcal{L}^{1,2}}\longrightarrow\mathcal{B}_{\bigtriangledown}$
is consistent with
$\mathcal{L}_{A}:\mathcal{L}^{2}(\mathbb{R})\longrightarrow\mathcal{L}^{2}(\mathbb{R})$.

###### Theorem 3.2

The linear canonical transform
$\mathcal{L}_{A}:\mathcal{B}_{\mathcal{L}^{1,2}}\longrightarrow\mathcal{B}_{\bigtriangledown}$
is a bijection.

###### Theorem 3.3

Let $F,G\in\mathcal{B}_{\mathcal{L}^{1,2}}$ then

  1. (a)

$\quad\mathcal{L}_{A}[F+\lambda
G]=\mathcal{L}_{A}(F)+\lambda\mathcal{L}_{A}(G)$, for any complex $\lambda$.

  2. (b)

$\quad\mathcal{L}_{A}[e^{ikt}F](u)=e^{\frac{-idk(2u-bk)}{2}}\mathcal{L}_{A}[F](u-bk)$,
for $k\in\mathbb{R}$.

  3. (c)

$\quad\mathcal{L}_{A}[F(t+\tau)](u)=e^{i(2u+a\tau)\frac{\tau}{2b}}\mathcal{L}_{A}[e^{\frac{-ia}{b}x\tau}F(x)](u)$.

  4. (d)

$\quad\mathcal{L}_{A}[F^{(2)}](u)=\bigg{[}\bigg{(}\frac{iu}{b}\bigg{)}^{2}+\frac{ia}{b}\bigg{]}\mathcal{L}_{A}[F(t)](u).$

###### Theorem 3.4

Let $F,G\in\mathcal{B}_{\mathcal{L}^{1,2}}$ then
$\mathcal{L}_{A}(F*^{A}G)=\mathcal{L}_{A}(F)\mathcal{L}_{A}(G)$.

###### Theorem 3.5

Let $\delta-\lim F_{n}=F$ for $F_{n},F\in\mathcal{B}_{\mathcal{L}^{1,2}}$ then
$\mathcal{L}_{A}(F_{n})\rightarrow\mathcal{L}_{A}(F)$ uniformly on each
compact set of $\mathbb{R}$.

Proof: Let $\\{\delta_{m}\\}$ be a delta sequence such that
$F_{n}*^{A}\delta_{m},F*^{A}\delta_{m}\in\mathcal{L}^{1,2}(\mathbb{R})$ for
all $n,m\in\mathbb{N}$ and $\|(F_{n}-F)*^{A}\delta_{m}\|_{2}\rightarrow 0$ as
$n\rightarrow\infty$ for each $m\in\mathbb{N}$. Let $M$ be a compact set in
$\mathbb{R}$ then $\mathcal{L}_{A}(\delta_{m})>0$ on $M$ for all most
$m\in\mathbb{N}$. Since $\mathcal{L}_{A}(\delta_{m})$ is a continuous function
and
$\mathcal{L}_{A}(F_{n})*^{A}\mathcal{L}_{A}(\delta_{m})-\mathcal{L}_{A}(F)*^{A}\mathcal{L}_{A}(\delta_{m})=((\mathcal{L}_{A}(F_{n})-\mathcal{L}_{A}(F))*^{A}\mathcal{L}_{A}(\delta_{m}))$,
implies
$\|(\mathcal{L}_{A}(F_{n})-\mathcal{L}_{A}(F))*^{A}\mathcal{L}_{A}(\delta_{m})\|_{2}\rightarrow
0$ as $n\rightarrow\infty$ for each $m\in\mathbb{N}$. Thus
$\mathcal{L}_{A}(F_{n})\rightarrow\mathcal{L}_{A}(F)$ uniformly on $M$.
$\hfill\blacksquare$

## References

  * [1] T. K. Boehme; The support of Mikusinski operators, Trans. Amer. Math. Soc., 176, 319-334, (1973).
  * [2] Deng Bing, Tao Ran and Wang Yue. Convolution theorems for the linear canonical transform and their applications. Sci. China Series F: Inf. Sci., 49(5), 592-603, (2006).
  * [3] Deyun Wei, Qiwen Ran and Yong Li; New convolution theorem for the linear canonical transform and its translation invariance property, Optik 123, 1478-1481, (2012).
  * [4] Pravinkumar V. Dole and S. K. Panchal, Linear canonical transform for Integrable Boehmians, Int. J. Pure Appl. Math., 116(1), 91-96, (2017).
  * [5] J.Mikusinski and P.Mikusinski; Quotients de suites et leurs applications dans l’analyse fonctionnelle, C.R. Acad. Sci. Paris Ser. I Math., 293, 463-464, (1981).
  * [6] P. Mikusinski; Convergence of Boehmianes, Japan. J. Math, 9(1), 169-179,(1983).
  * [7] P. Mikusinski; Fourier transform for integrable Boehmians, Rocky Moun-tain J. Math., 17(3), 577-582,(1987).
  * [8] P. Mikusinski; Boehmians and generalized functions, Acta. Math. Hungarica, 51, 159-179, (1988).
  * [9] P. Mikusinski; Transform of Boehmians, Different Aspects of Differentiability Dissertationes Mathematicae, 340, 201-206,(1995).
  * [10] M. Moshinsky and C. Quesne; Linear canonical transformations and their unitary representations, J. Math. Phys, 12(8), 1772-1783,(1971).
  * [11] Dennis Nemzer; Laplace transforms on a class of Boehmians, Bull. Austral. Math. Soc., 46, 347-352, (1992).
  * [12] R. S. Pathak, _A Course in Distributional Theory and Applications_ , Narosa Publication House, New Delhi 2001.
  * [13] A. Singh and P. K. Banergi; Fractional integrals of fractional Fourier transform for integrable Boehmians, Proceedings of the National Academy of Sciences, India Section A: Physical Sciences, 2017.
  * [14] Walter Rudin; Real and Complex Analysis, Third Edition, McGraw-Hill, New York, 1987.
  * [15] A. I. Zayed; Fractional Fourier transform of generalized functions, Integ. Trans. Spl. Funct.,7, 299-312, (1998).


# Multiple importance sampling for stochastic gradient estimation

Corentin Salaün  Xingchang Huang  Iliyan Georgiev  Niloy J. Mitra  Gurprit
Singh

###### Abstract

We introduce a theoretical and practical framework for efficient importance
sampling of mini-batch samples for gradient estimation from single and
multiple probability distributions. To handle noisy gradients, our framework
dynamically evolves the importance distribution during training by utilizing a
self-adaptive metric. Our framework combines multiple, diverse sampling
distributions, each tailored to specific parameter gradients. This approach
facilitates the importance sampling of _vector-valued_ gradient estimation.
Rather than naively combining multiple distributions, our framework involves
optimally weighting data contribution across multiple distributions. This
adapted combination of multiple importance yields superior gradient estimates,
leading to faster training convergence. We demonstrate the effectiveness of
our approach through empirical evaluations across a range of optimization
tasks like classification and regression on both image and point cloud
datasets.

Machine Learning, ICML

##  1 Introduction

Stochastic gradient descent (SGD), in tandem with gradient backpropagation, is
fundamental in optimizing complex neural networks. This iterative optimization
process relies on the efficient estimation of gradients to update model
parameters and minimize the optimization objective. A significant challenge in
methods based on SGD lies in the influence of stochasticity on gradient
estimation, impacting both the quality of the estimates and convergence speed.
This stochasticity introduces errors in the form of noise, and addressing and
minimizing such noise in gradient estimation continues to be an active area of
research.

Various approaches have been introduced to reduce gradient estimation noise,
including data diversification (Zhang et al., 2019; Faghri et al., 2020; Ren
et al., 2019), adaptive mini-batch sizes (Balles et al., 2017; Alfarra et al.,
2021), momentum-based estimation (Rumelhart et al., 1986; Kingma & Ba, 2014),
and adaptive sampling strategies (Santiago et al., 2021). These methods
collectively expedite the optimization by improving the gradient-estimation
accuracy.

Another well-established technique for noise reduction in estimation is
importance sampling (IS) (Loshchilov & Hutter, 2015; Katharopoulos & Fleuret,
2017, 2018), which involves the non-uniform selection of data samples for
mini-batch construction. Data samples that contribute more significantly to
gradient estimation are selected more often. This allows computational
resources to focus on the most critical data for the optimization task.
However, these algorithms are quite inefficient and add significant overhead
to the training process. Another limitation of importance sampling, in
general, lies in determining the best sampling distribution to achieve maximal
improvement, often necessitating a quality trade-off due to the simultaneous
estimation of numerous parameters.

We propose an efficient importance sampling algorithm that does not require
resampling, in contrast to (Katharopoulos & Fleuret, 2018). Our importance
function dynamically evolves during training, utilizing a self-adaptive metric
to effectively manage initial noisy gradients. Further, unlike existing IS
methods in machine learning where importance distributions assume scalar-
valued gradients, we propose a multiple importance sampling (MIS) strategy to
manage _vector-valued_ gradient estimation. We propose the simultaneous use of
multiple sampling strategies combined with a weighting approach following the
principles of MIS theory, well studied in the rendering literature in computer
graphics (Veach, 1997). Rather than naively combining multiple distributions,
our proposal involves estimating importance weights w.r.t. data samples across
multiple distributions by leveraging the theory of optimal MIS (OMIS)
(Kondapaneni et al., 2019). This optimization process yields superior gradient
estimates, leading to faster training convergence. In summary, we make the
following contributions:

  * •

We develop an efficient IS algorithm with a self-adaptive metric for
importance sampling.

  * •

We introduce an MIS estimator for gradient estimation to importance sample
vector-valued gradients for gradient estimation.

  * •

We present a practical approach to compute the optimal weights for multiple
sampling strategies to maximize the quality of vector-valued gradient
estimation.

  * •

We demonstrated the effectiveness of our approach on various machine learning
tasks.

##  2 Related work

|  |  | Output 1 | Output 2 | Output 3
---|---|---|---|---|---
|  |  |  |  |
|  |  |  |
(a) Network diagram | (b) Ground-truth | (c) Output-layer | (d) Norms of individual output nodes
| classification | gradient norm |  |  |
Figure 1:  We visualize different importance sampling distributions for a
simple classification task. We propose to use the output layer gradients for
importance sampling, as shown in the network diagram (a). For a given ground-
truth classification (top) and training dataset (bottom) shown in (b), it is
possible to importance sample from the $L_{2}$ norm of the output-layer
gradients (c) or from three different sampling distributions derived from the
gradient norms of individual output nodes (d). The bottom row shows sample
weights from each distribution.

### Importance sampling for gradient estimation.

Importance sampling (IS) (Kahn, 1950; Kahn & Marshall, 1953; Owen & Zhou,
2000) has emerged as a powerful technique in high energy physics, Bayesian
inference, rare event simulation for finance and insurance, and rendering in
computer graphics. In the past few years, IS has also been applied in machine
learning to improve the accuracy of gradient estimation and enhance the
overall performance of learning algorithms (Zhao & Zhang, 2015).

By strategically sampling data points from a non-uniform distribution, IS
effectively focuses training resources on the most informative and impactful
data, leading to more accurate gradient estimates. Bordes et al. (2005)
developed an online algorithm (LASVM) that uses importance sampling to train
kernelized support vector machines. Loshchilov & Hutter (2015) suggested
employing data rankings based on their respective loss values. This ranking is
then employed to create an importance sampling strategy that assigns greater
importance to data with higher loss values. Katharopoulos & Fleuret (2017)
proposed importance sampling the loss function. Subsequently, Katharopoulos &
Fleuret (2018) introduced an upper bound to the gradient norm that can be
employed as an importance function. Their algorithm involves resampling and
computing gradients with respect to the final layer. Despite the importance
function demonstrating improvement over uniform sampling, their algorithm
exhibits significant inefficiency.

Appendix B summarizes the theory behind (multiple) importance sampling. It
also states the optimal MIS estimator and how to compute it.

### Multiple importance sampling.

The concept of Multiple Importance Sampling (MIS) emerged as a robust and
efficient technique for integrating multiple sampling strategies (Owen & Zhou,
2000). Its core principle lies in assigning weights to various importance
sampling estimation, allowing each data sample to utilize the most appropriate
strategy. Veach (1997) introduced this concept of MIS to rendering in computer
graphics and proposed the widely adopted _balance heuristic_ for importance
(weight) allocation. The balance heuristic determines weights based on a data
sample’s relative importance across all sampling approaches, effectively
mitigating the influence of outliers with low probability densities. While MIS
is straightforward to implement and independent of the specific function,
Variance-Aware MIS (Grittmann et al., 2019) advanced the concept by using
variance estimates from each sampling technique for further error reduction.
Moreover, Optimal MIS (Kondapaneni et al., 2019) derived optimal sampling
weights that minimize MIS estimator variance. Notably, these weights depend
not only on probability density but also on the function values of the
samples.

##  3 Problem statement

The primary goal of machine-learning optimization is to find the optimal
parameters $\theta$ for a given model function $m(x,\theta)$ by minimizing a
loss function ${\mathcal{L}}$ over a dataset ${\Omega}$:

$\displaystyle\theta^{*}=\underset{\theta}{\mathrm{argmin}}\,\underbrace{\int_{{\Omega}}{\mathcal{L}}(m(x_{i},\theta),y)\,\mathrm{d}x.}_{L_{\theta}}$
(1)

The loss function ${\mathcal{L}}$ quantifies the dissimilarity between the
model predictions $m(x,\theta)$ and observed data $y$. The factor in front of
the integral normalizes the overall loss $L_{\theta}$ with respect to the
dataset size. In the common case of a discrete dataset, the integral becomes a
sum.

In practice, the total loss is minimized via iterative gradient descent. In
each iteration $t$, the gradient $\nabla L_{\theta_{t}}$ of the loss with
respect to the current model parameters $\theta_{t}$ is computed, and the
parameters are updated as

$\theta_{t+1}=\theta_{t}-\lambda\underbrace{\int_{{\Omega}}\nabla{\mathcal{L}}(m(x,\theta),y)\,\mathrm{d}x}_{\nabla
L_{\theta_{t}}},$ (2)

where $\lambda>0$ is the learning rate.

### Monte Carlo gradient estimator.

In practice, the parameter gradient is estimated from a small batch
$\\{x_{i}\\}_{i=1}^{B}$ of randomly selected data points:

$\langle\nabla
L_{\theta}\rangle=\sum_{i=1}^{B}\frac{\nabla{\mathcal{L}}(m(x_{i},\theta),y_{i})}{{B}p(x_{i})}\approx\nabla
L_{\theta},\quad x_{i}\sim p.$ (3)

The data points are sampled from a probability density function (pdf) $p$ or
probability mass function in discreet cases. The mini-batch gradient descent
substitutes the true gradient $\nabla L_{\theta_{t}}$ with an estimate
$\langle\nabla L_{\theta_{t}}\rangle$ in Equation 2 to update the model
parameters in each iteration.

We want to estimate $\nabla L_{\theta_{t}}$ accurately and also efficiently,
since the gradient-descent iteration (2) may require many thousands of
iterations until the parameters converge. These goals can be achieved by
performing the optimization in small batches whose samples are chosen
according to a carefully designed distribution $p$. For a simple
classification problem, Figure 1c shows an example importance sampling
distribution derived from the output layer of the model. In Figure 1d we
derive multiple distributions from the individual output nodes. Below we
develop theory and practical algorithms for importance sampling using a single
distribution (Section 4) and for combining multiple distributions to further
improve gradient estimation (Section 5).

##  4 Mini-batch importance sampling

Mini-batch gradient estimation (3) notoriously suffers from Monte Carlo noise,
which can make the parameter-optimization trajectory erratic and convergence
slow. That noise comes from the often vastly different contributions of
different samples $x_{i}$ to that estimate.

Typically, the selection of samples that go into a mini-batch is done with
uniform probability $p(x_{i})=1/{|\Omega|}$. Importance sampling is a
technique for using a non-uniform pdf to strategically pick samples
proportionally on their contribution to the gradient, to reduce estimation
variance.

### Practical algorithm.

We propose an importance sampling algorithm for mini-batch gradient descent,
outlined in Algorithm 1. Similarly to Schaul et al. (2015), we use an
importance function that relies on readily available quantities for each data
point, introducing only negligible memory and computational overhead over
classical uniform mini-batching. We store a set of persistent _un-normalized
importance_ scalars $q=\\{q_{i}\\}_{i=1}^{|\Omega|}$ that are updated
continuously during the optimization.

The first epoch is a standard SGD one, during which we additionally compute
the initial importance of each data point (line 3). In each subsequent epoch,
at each mini-batch optimization step $t$ we normalize the importance values to
a valid distribution $p$ (line 6). We then choose ${B}$ data samples (with
replacement) according to $p$ (line 7). The loss ${\mathcal{L}}$ is evaluated
for each selected data sample (line 8), and backpropagated to compute the loss
gradient (line 9). The per-sample importance is used in the gradient
estimation (line 10) to normalize the contribution. In practice lines 9-10 can
be done simultaneously by backpropagating a weighted loss
${\mathcal{L}}(x)\cdot(\nicefrac{{1}}{{(p(x)\cdot B)}})^{T}$. Finally, the
network parameters are updated using the estimated gradient (line 11). On line
12, we update the importance of the samples in the mini-batch; we describe our
choice of importance function below. The blending parameter $\gamma$ ensures
stability of the persistent importance as discussed in Appendix E. At the end
of each epoch (line 13), we add a small value to the un-normalized weights of
all data to ensure that every data point will be eventually evaluated, even if
its importance is deemed low by the importance metric.

Algorithm 1 Mini-batch importance sampling for SGD.

1:$\theta\leftarrow$ random parameter initialization

2:$B\leftarrow$ mini-batch size, $N=|\Omega|$ $\leftarrow$ Dataset size

3:$q,\theta\leftarrow\text{Initialize}(x,y,\Omega,\theta,B)$ $\leftarrow$
Algorithm 4

4:until convergence do $\leftarrow$ Loop over epochs

5: for $t\leftarrow 1$ to $N/B$ do $\leftarrow$ Loop over mini-batches

6: $p\leftarrow q/$sum$(q)$ $\leftarrow$ Normalize importance to pdf

7: $x,y\leftarrow{B}$ data samples $\\{x_{i},y_{i}\\}_{i=1}^{B}\propto p$

8: ${\mathcal{L}}(x)\leftarrow{\mathcal{L}}(m(x,\theta),y)$

9: $\nabla{\mathcal{L}}(x)\leftarrow$ Backpropagate$({\mathcal{L}}(x))$

10: $\langle\nabla
L_{\theta}\rangle\leftarrow(\nabla{\mathcal{L}}(x)\cdot(\nicefrac{{1}}{{p(x)}})^{T})/B$
$\leftarrow$ Equation 3

11: $\theta\leftarrow\theta-\eta\,\langle\nabla L_{\theta}\rangle$
$\leftarrow$ SGD step

12: $q(x)\leftarrow\gamma\cdot
q(x)+(1-\gamma)\cdot\left\|\frac{\partial{\mathcal{L}}(x)}{\partial
m(x,\theta)}\right\|$

13: $q\leftarrow q+\epsilon$ $\rcurvearrowne$ Accumulate importance

14:return $\theta$

It is important to note that the first epoch is done without importance
sampling to initialize each sample importance. This does not add overhead as
it is equivalent to a classical epoch running over all data samples. While
similar schemes have been proposed in the past (Loshchilov & Hutter, 2015),
they often rely on a multitude of hyperparameters, making their practical
implementation challenging. This has led to the development of alternative
methods like re-sampling (Katharopoulos & Fleuret, 2018; Dong et al., 2021;
Zhang et al., 2023). Tracking importance across batches and epochs minimizes
the computational overhead, further enhancing the efficiency and practicality
of the approach.

### Importance function.

In combination with the presented algorithm, we propose an importance function
that is efficient to evaluate. While the gradient $L_{2}$ norm has been shown
to be optimal (Zhao & Zhang, 2015; Needell et al., 2014; Wang et al., 2017;
Alain et al., 2015), calculating it can be computationally expensive as it
requires full backpropagation for every data point. To this end, we compute
the gradient norm only for a subset of the parameters, specifically the output
nodes of the network: $q(x)=\left\|\frac{\partial\mathcal{L}(x)}{\partial
m(x,\theta)}\right\|$. This choice is based on an upper bound of the gradient
norm, using the chain rule and the Cauchy–Schwarz inequality (Katharopoulos &
Fleuret, 2018):

$\displaystyle\\!\left\|\frac{\partial\mathcal{L}(x_{i})}{\partial\theta}\right\|$
$\displaystyle=\left\|\frac{\partial\mathcal{L}(x)}{\partial
m(x,\theta)}\cdot\frac{\partial m(x,\theta)}{\partial\theta}\right\|\leq$ (4)
$\displaystyle\left\|\frac{\partial\mathcal{L}(x)}{\partial
m(x,\theta)}\right\|\cdot\left\|\frac{\partial
m(x,\theta)}{\partial\theta}\right\|\leq\underbrace{\left\|\frac{\partial\mathcal{L}(x)}{\partial
m(x,\theta)}\right\|}_{q(x)}\cdot\,C\,,$

where $C$ is the Lipschitz constant of the parameters gradient. That is, our
importance function is a bound of the gradient magnitude based on the output-
layer gradient norm.

We tested the relationship between four different importance distributions:
uniform, our proposed importance function,

the loss function as importance (Katharopoulos & Fleuret, 2017), and the work
by Katharopoulos & Fleuret (2018) using an other gradient norm bound. The
inline figure plots the $L_{2}$ difference between these importance
distributions and the ground-truth gradient-norm distribution across epochs
for an MNIST classification task. It shows that Our IS distribution has the
smallest difference, i.e., it achieves high accuracy while requiring only a
small part of the gradient.

For some specific task when the output layer has predictable shape, it is
possible to derive a closed form definition of the proposed importance metric.
Appendix D derives the close form importance for classification task using
cross entropy loss.

Note that any importance heuristic can be used on line 12 of Algorithm 1, such
as the gradient norm (Zhao & Zhang, 2015; Needell et al., 2014; Wang et al.,
2017; Alain et al., 2015), the loss (Loshchilov & Hutter, 2015; Katharopoulos
& Fleuret, 2017; Dong et al., 2021), or more advanced importance
(Katharopoulos & Fleuret, 2018). For efficiency, our importance function
reuses the forward-pass computations from line 8, updating $q$ only for the
current mini-batch samples.

##  5 Multiple importance sampling

The parameter gradient $\nabla L_{\theta}$ is vector with dimension equal to
the number of model parameters. The individual parameter derivatives vary
uniquely across the data points, and estimation using a single distribution
(Section 4) inevitably requires making a trade-off, e.g., only importance
sampling the overall gradient magnitude. Truly minimizing the estimation error
requires estimating each derivative using a separate importance sampling
distribution tailored to its variation. However, there are two practical
issues with this approach: First, it would necessitate sampling from all of
these distributions, requiring “mini-batches” of size equal at least to the
number of parameters. Second, it would lead to significant computation waste,
since backpropagation computes all parameter derivatives but only one of them
would be used per data sample. To address this issue, we propose using a small
number of distributions, each tailored to the variation of a parameter subset,
and combining _all_ computed derivatives into a low-variance estimator, using
multiple importance sampling theory. As an example, Figure 1d shows three
sampling distributions for a simple classification task, based on the
derivatives of the network’s output nodes, following the boundary of each
class.

### MIS gradient estimator.

Combining multiple sampling distributions into a single robust estimator has
been well studied in the Monte Carlo rendering literature. The best known
method is _multiple importance sampling_ (MIS) (Veach, 1997). In our case of
gradient estimation, the MIS estimator takes for form

$\langle\nabla
L_{\theta}\rangle_{\mathrm{MIS}}=\sum_{j=1}^{J}\sum_{i=1}^{n_{j}}w_{j}(x_{ij})\frac{\nabla{\mathcal{L}}(m(x_{ij},\theta),y_{ij})}{n_{j}p_{j}(x_{ij})},$
(5)

where $J$ is the number of sampling distributions, $n_{j}$ the number of
samples from distribution $j$, and $x_{ij}$ the $i$th sample from the $j$th
distribution. Each sample is modulated by a weight $w_{j}(x_{ij})$; the
estimator is unbiased as long as $\sum_{j=1}^{J}w_{j}(x)=1$ for every data
point $x$ in the dataset.

### Optimal weighting.

Various MIS weighting functions $w_{j}$ have been proposed in literature, the
most universally used one being the balance heuristic (Veach, 1997). In this
work we use the recently derived optimal weighting scheme (Kondapaneni et al.,
2019) which minimizes the estimation variance for a given set of sampling
distributions $p_{j}$:

$w_{j}(x)=\alpha_{j}\frac{p_{j}(x)}{\nabla{\mathcal{L}}(m(x,\theta),y)}\;+\\\
\frac{n_{j}p_{j}(x)}{\sum_{k=1}^{J}n_{k}p_{k}(x)}\Bigg{(}1-\frac{\sum_{k=1}^{J}\alpha_{k}p_{k}(x)}{\nabla{\mathcal{L}}(m(x,\theta),y)}\Bigg{)}.$
(6)

Here, $\boldsymbol{\alpha}=[\alpha_{1},\ldots,\alpha_{J}]$ is the solution to
the linear system

$\small\boldsymbol{A}\boldsymbol{\alpha}=\boldsymbol{b}\text{, with
}\begin{dcases}a_{j,k}=\int_{{\Omega}}\frac{p_{j}(x)p_{k}(x)}{\sum_{i}^{J}n_{i}p_{i}(x)}d(x,y),\\\
b_{j}=\int_{{\Omega}}\frac{p_{j}(x)\nabla{\mathcal{L}}(m(x,\theta),y)}{\sum_{i}^{J}n_{i}p_{i}(x)}d(x,y),\end{dcases}$
(7)

where $a_{j,k}$ and $b_{j}$ are the elements of the matrix
$\boldsymbol{A}\in\mathbb{R}^{J\times J}$ and vector
$\boldsymbol{b}\in\mathbb{R}^{J}$ respectively.

Instead of explicitly computing the optimal weights in Equation 6 using
Equation 7 and plugging them into the MIS estimator (5), we can use a shortcut
evaluation that yields the same result (Kondapaneni et al., 2019):

$\langle\nabla L_{\theta}\rangle_{\mathrm{OMIS}}=\sum_{j=1}^{J}\alpha_{j}.$
(8)

In Appendix B we provide an overview of MIS and the aforementioned weighting
schemes. Importantly for our case, the widely adopted balance heuristic does
not bring practical advantage over single-distribution importance sampling
(Section 4) as it is equivalent to sampling from a mixture of the given
distributions; we can easily sample from this mixture by explicitly averaging
the distributions into a single one. In contrast, the optimal weights are
different for each gradient dimension as they depend on the gradient value.

Algorithm 2 Optimal multiple importance sampling SGD.

1:$\theta\leftarrow$ random parameter initialization

2:$B\leftarrow$ mini-batch size, $J\leftarrow$ number of pdf

3:$N=|\Omega|\leftarrow$ dataset size

4:$n_{j}\leftarrow$ sample count per technique, for $j\in\\{1,..J\\}$

5:$\boldsymbol{q},\theta\leftarrow\text{InitializeMIS}(x,y,\Omega,\theta,B)$
$\leftarrow$ Algorithm 5

6:$\langle\boldsymbol{A}\rangle\leftarrow 0^{J\times
J},\langle\boldsymbol{b}\rangle\leftarrow 0^{J}$ $\leftarrow$ OMIS linear
system

7:until convergence do $\leftarrow$ Loop over epochs

8: for $t\leftarrow 1$ to $N/B$ do $\leftarrow$ Loop over mini-batches

9:
$\langle\boldsymbol{A}\rangle\leftarrow\beta\langle\boldsymbol{A}\rangle,\langle\boldsymbol{b}\rangle\leftarrow\beta\langle\boldsymbol{b}\rangle$

10: for $j\leftarrow 1$ to $J$ do $\leftarrow$ Loop over distributions

11: $p_{j}\leftarrow q_{j}/\text{sum}(q_{j})$

12: $x,y\leftarrow{B}$ data samples $\\{x_{i},y_{i}\\}_{i=1}^{n_{j}}\propto
p_{j}$

13: ${\mathcal{L}}(x)\leftarrow{\mathcal{L}}(m(x,\theta),y)$

14: $\nabla{\mathcal{L}}(x)\leftarrow$ Backpropagate$({\mathcal{L}}(x))$

15: $S(x)\leftarrow\sum_{k=1}^{J}n_{k}p_{k}(x)$

16:
$\boldsymbol{W}\leftarrow\nicefrac{{n_{i}p_{i}(x)}}{{\sum_{k=1}^{J}n_{k}p_{k}(x)}}$
$\rcurvearrowsw$ Momentum estim.

17:
$\langle\boldsymbol{A}\rangle\leftarrow\langle\boldsymbol{A}\rangle+(1-\beta)\sum_{i=1}^{n_{j}}\boldsymbol{W}_{i}\boldsymbol{W}_{i}^{T}$

18:
$\langle\boldsymbol{b}\rangle\leftarrow\langle\boldsymbol{b}\rangle+(1-\beta)\sum_{i=1}^{n_{j}}\nabla{\mathcal{L}}(x_{i})\boldsymbol{W}_{i}/S(x_{i})$

19:
$\boldsymbol{q}(x)\leftarrow\gamma\boldsymbol{q}(x)+(1-\gamma)\frac{\partial\mathcal{L}(x)}{\partial
m(x,\theta)}$

20:
$\langle\boldsymbol{\alpha}\rangle\leftarrow\langle\boldsymbol{A}\rangle^{-1}\langle\boldsymbol{b}\rangle$

21: $\langle\nabla
L_{\theta}\rangle_{\mathrm{OMIS}}\leftarrow\sum_{j=1}^{J}\langle\alpha_{j}\rangle$

22: $\theta\leftarrow\theta-\eta\,\langle\nabla
L_{\theta}\rangle_{\mathrm{OMIS}}$ $\leftarrow$ SGD step

23:return $\theta$

### Practical algorithm.

Implementing the optimal-MIS estimator (8) amounts to drawing $n_{j}$ samples
from each distribution, computing $\boldsymbol{\alpha}$ for each dimension of
the gradient and summing its elements. The integrals in $\boldsymbol{A}$ and
$\boldsymbol{b}$ (sums in the discrete-dataset case) can be estimated as
$\langle\boldsymbol{A}\rangle$ and $\langle\boldsymbol{b}\rangle$ from the
drawn samples, yielding the estimate
$\langle\boldsymbol{\alpha}\rangle=\langle\boldsymbol{A}\rangle^{-1}\langle\boldsymbol{b}\rangle$.

Algorithm 2 shows a complete gradient-descent algorithm. The main differences
with Algorithm 1 are the use of multiple importance distributions
$\boldsymbol{q}=\\{q_{j}\\}_{j=1}^{J}$ (line 5) and the linear system used to
compute the OMIS estimator (line 6). This linear system is updated (lines
15-18) using the mini-batch samples and solved to obtain the gradient
estimation (line 21). Since the matrix $\langle\boldsymbol{A}\rangle$ is
independent of the gradient estimation (see Equation 7), its inversion can be
shared across all parameter estimates.

Figure 2:  Convergence comparison of polynomial regression of order 6 using
different method. Exact gradient show a gradient descent as baseline and
classical SGD. For our method, we compare importance sampling and OMIS using
$n=2$ or $4$ importance distributions. Balance heuristic MIS is also visible.
Our method using OMIS achieve same convergence as exact gradient.

### Momentum-based linear-system estimation.

If the matrix estimate $\langle\boldsymbol{A}\rangle$ is inaccurate, its
inversion can be unstable and yield a poor gradient estimate. The simplest way
to tackle this problem is to use a large number of samples per distribution,
which produces a accurate estimates of both $\boldsymbol{A}$ and
$\boldsymbol{b}$ and thus a stable solution to the linear system. However,
this approach is computationally expensive. Instead, we keep the sample counts
low and reuse the estimates from previous mini-batches via momentum-based
accumulation, shown in lines 17–18, where $\beta$ is the parameter controlling
the momentum; we use $\beta=0.7$. This accumulation provides stability, yields
an estimate of the momentum gradient (Rumelhart et al., 1986), and allows us
to use 1–4 samples per distribution in a mini-batch.

### Importance functions.

To define our importance distributions, we expand on the approach from Section
4. Instead of taking the norm of the entire output layer of the model, we take
the different gradients separately as
$\boldsymbol{q}(x)=\frac{\partial\mathcal{L}(x)}{\partial m(x,\theta)}$ (see
Figure 1d). Similarly to Algorithm 1, we apply momentum-based accumulation of
the per-data importance (line 19 in Algorithm 2). If the output layer has more
nodes than the desired number $J$ of distributions, we select a subset of the
nodes. Many other ways exist to derive the distributions, e.g., clustering the
nodes into $J$ groups and taking the norm of each; we leave such exploration
for future work.

##  6 Experiments

### Implementation details.

We evaluate our importance sampling (IS) and optimal multiple importance
sampling (OMIS) methods on a set of classification and regression tasks with
different data modalities (images, point clouds). We compare them to classical
SGD (which draws mini-batch samples uniformly without replacement), DLIS
(Katharopoulos & Fleuret, 2018), and LOW (Santiago et al., 2021). DLIS uses a
resampling scheme that samples an initial, larger mini-batch uniformly and
then selects a fraction of them for backpropagation and a gradient step. This
resampling is based on an importance sampling metric computed by running a
forward pass for each initial sample. LOW applies adaptive weighting to
uniformly selected mini-batch samples to give importance to data with high
loss. All reported metrics are computed on data unseen during training, with
the exception of the regression tasks.

All experiments are conducted on a single NVIDIA Tesla A40 graphics card.
Details about the optimization setup of each experiment can be found in
Appendix A.

### Convex problem.

We performed a basic convergence analysis of IS and OMIS on a convex
polynomial-regression problem. Figure 2 compares classical SGD, our IS, and
three MIS techniques: balance heuristic (Veach, 1997) and our OMIS using two
and four importance distributions. The exact gradient serves as a reference
point for optimal convergence. Balance-heuristic MIS exhibits similar
convergence to IS. This can be attributed to the weights depending solely on
the relative importance distributions, disregarding differences in individual
parameter derivatives. This underscores the unsuitability of the balance
heuristic as a weighting method for vector-valued estimation. Both our OMIS
variants achieve convergence similar to that of the exact gradient. The four-
distribution variant achieves the same quality as the exact gradient using
only 32 data samples per mini-batch. This shows the potential of OMIS to
achieve low error in gradient estimation even at low mini-batch sizes.

Figure 3:  Classification error convergence for MNIST classification for
various methods. Both Katharopoulos & Fleuret (2018) (DLIS) and resampling SGD
approach. In comparison, our three method use the presented algorithm without
resampling. It is visible that while DLIS perform similarly to our IS at equal
epoch, the overhead of the method makes ours noticeably better at equal time
for of IS and OMIS.

### Classification.

In Figure 3, we compare our algorithms to the DLIS resampling algorithm of
Katharopoulos & Fleuret (2018) on MNIST classification. Our IS performs
slightly better than DLIS, and our OMIS does best. The differences between our
methods and the rest are more pronounced when comparing equal-time
performance. DLIS has a higher computational cost as it involves running a
forward pass on a large mini-batch to compute resampling probabilities. Our
OMIS requires access to the gradient of each mini-batch sample; obtaining
these gradients in our current implementation is inefficient due to technical
limitations in the optimization framework we use (PyTorch). Nevertheless, the
method manages to make up for this overhead with a higher-quality gradient
estimate. In Figure 3 we compare classification error; loss-convergence plots
are shown in Appendix F (Figure 8).

Figure 4:  On CIFAR-100 classification dataset, instead of comparing the DLIS
resampling algorithm , we use DLIS importance metric in our Algorithm 1. We
display zoom-in of the end of the curves to highlight the differences. At
equal epochs (left), our methods (Our IS & Our AS) show improvements compared
to LOW (Santiago et al., 2021) and DLIS weights. At equal time (right), LOW
and the DLIS weights takes longer to converge.Overall our approach shows
faster convergence with lower importance computation.  Figure 5:  Comparisons
on CIFAR-10 using Vision Transformer (ViT) (Dosovitskiy et al., 2020). The
results show our importance sampling scheme (Our IS) can improve over
classical SGD, LOW (Santiago et al., 2021) and DLIS (Katharopoulos & Fleuret,
2018) on modern transformer architecture.

In Figure 4, we compare our IS against using the DLIS importance function in
Algorithm 1 and LOW (Santiago et al., 2021) on CIFAR-100 classification. At
equal number of epochs, the difference between the methods is small (see
close-up view). Our IS achieves similar classification accuracy as LOW and
outperforms the DLIS variant. At equal time the difference is more important
as our method has lower computational cost. This experiment shows that our
importance function achieves better performance than that of DLIS within the
same optimization algorithm.

Figure 5 shows a similar experiment on CIFAR-10 using a vision transformer
(Dosovitskiy et al., 2020). Our IS method achieves consistent improvement over
the state of the art. The worse convergence of (original, resampling-based)
DLIS can be attributed to its resampling tending to exclude some training data
with very low importance, which can cause overfitting.

Figure 6:  Comparison of our two methods (Our IS, Our OMIS) on point-cloud classification using PointNet (Qi et al., 2017) architecture. Our OMIS achieves lower classification error at equal epochs, though it introduces computation overhead as shown at equal-time comparisons. At equal time, our method using importance sampling achieves the best performance.  |  |  |  |  |
---|---|---|---|---|---
|  |  |  |
Reference | Uniform | DLIS | Our IS | Our OMIS
Figure 7:  Comparison at equal step for image 2D regression. Left side show
the convergence plot while the right display the result regression and a
close-up view. Our method using MIS achieves the lower error on this problem
while IS and DLIS perform similarly. On the images it is visible that our OMIS
recover the finest details of the fur and whiskers.

Figure 6 shows point-cloud classification, where our IS is comparable to
classical SGD and our OMIS outperforms other methods in terms of
classification error at equal epochs. Equal-time comparison demonstrates that
our IS is as efficient an SGD in complexes cases where importance sampling
does not improve convergence. DLIS and our OMIS both suffer from computational
overhead.

We also perform an ablation study for linear-system momentum in Algorithm 2.
We apply same momentum on the gradient for classical SGD, DLIS and our IS.
Appendix F (Figure 9) shows this comparison. Our OMIS still outperforms other
methods for this task at equal steps.

### Regression.

Figure 7 shows results on image regression, comparing classical SGD, DLIS, and
our IS and OMIS. Classical SGD yields a blurry image, as seen in the zoom-ins.
DLIS and our IS methods achieves similar results, with increased whisker
sharpness but still blurry fur, though ours has slightly lower loss and is
computationally faster, as discussed above. Our OMIS employs three sampling
distributions based on the network’s outputs which represent the red, green
and blue image channels. This method achieves the lowest error and highest
image fidelity, as seen in the zoom-in.

##  7 Limitations and future work

We have showcased the effectiveness of importance sampling and optimal
multiple importance sampling (OMIS) in machine-learning optimization, leading
to a reduction in gradient-estimation error. Our current OMIS implementation
incurs some overhead as it requires access to individual mini-batch sample
gradients. Modern optimization frameworks can efficiently compute those
gradients in parallel but only return their average. This is the main
computational bottleneck in the method. The overhead of the linear system
computation is negligible; we have tested using up to 10 distributions.

Our current OMIS implementation is limited to sequential models; hence its
absence from our ViT experiment in Figure 5. However, there is no inherent
limitation that would prevent its use with such more complex architectures. We
anticipate that similar improvements could be achieved, but defer the
exploration of this extension to future work.

In all our experiments we allocate the same sampling budget to each
distribution. Non-uniform sample distribution could potentially further reduce
estimation variance, especially if it can be dynamically adjusted during the
optimization process.

Recent work from Santiago et al. (2021) has explored a variant of importance
sampling that forgoes sample-contribution normalization, i.e., the division by
the probability $p(x)$ in Equation 3 (and on line 10 of Algorithm 1). This
heuristic approach lacks proof of convergence but can achieve practical
improvement over importance sampling in some cases. We include a such variant
of our IS method in Appendix F.

##  8 Conclusion

This work proposes a novel approach to improve gradient-descent optimization
through efficient data importance sampling. We present a method incorporates a
gradient-based importance metric that evolves during training. It boasts
minimal computational overhead while effectively exploiting the gradient of
the network output. Furthermore, we introduce the use of (optimal) multiple
importance sampling for vector-valued, gradient estimation. Empirical
evaluation on typical machine learning tasks demonstrates the tangible
benefits of combining several importance distributions in achieving faster
convergence.

## Ethic and impact

This paper presents work that aims to advance the field of Machine Learning.
There are many potential societal consequences of our work, none of which we
feel must be specifically highlighted here.

## References

  * Alain et al. (2015) Alain, G., Lamb, A., Sankar, C., Courville, A., and Bengio, Y.  Variance reduction in sgd by distributed importance sampling.  _arXiv preprint arXiv:1511.06481_ , 2015.
  * Alfarra et al. (2021) Alfarra, M., Hanzely, S., Albasyoni, A., Ghanem, B., and Richtarik, P.  Adaptive learning of the optimal batch size of sgd, 2021.
  * Balles et al. (2017) Balles, L., Romero, J., and Hennig, P.  Coupling adaptive batch sizes with learning rates.  In _Proceedings of the 33rd Conference on Uncertainty in Artificial Intelligence (UAI)_ , pp. ID 141, August 2017.  URL http://auai.org/uai2017/proceedings/papers/141.pdf.
  * Bordes et al. (2005) Bordes, A., Ertekin, S., Weston, J., and Bottou, L.  Fast kernel classifiers with online and active learning.  _Journal of Machine Learning Research_ , 6(54):1579–1619, 2005.  URL http://jmlr.org/papers/v6/bordes05a.html.
  * Deng (2012) Deng, L.  The mnist database of handwritten digit images for machine learning research [best of the web].  _IEEE signal processing magazine_ , 29(6):141–142, 2012.
  * Dong et al. (2021) Dong, C., Jin, X., Gao, W., Wang, Y., Zhang, H., Wu, X., Yang, J., and Liu, X.  One backward from ten forward, subsampling for large-scale deep learning.  _arXiv preprint arXiv:2104.13114_ , 2021.
  * Dosovitskiy et al. (2020) Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv preprint arXiv:2010.11929_ , 2020.
  * Faghri et al. (2020) Faghri, F., Duvenaud, D., Fleet, D. J., and Ba, J.  A study of gradient variance in deep learning.  _arXiv preprint arXiv:2007.04532_ , 2020.
  * Grittmann et al. (2019) Grittmann, P., Georgiev, I., Slusallek, P., and Křivánek, J.  Variance-aware multiple importance sampling.  _ACM Trans. Graph._ , 38(6), nov 2019.  ISSN 0730-0301.  doi: 10.1145/3355089.3356515.  URL https://doi.org/10.1145/3355089.3356515.
  * He et al. (2016) He, K., Zhang, X., Ren, S., and Sun, J.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770–778, 2016.
  * Kahn (1950) Kahn, H.  Random sampling (monte carlo) techniques in neutron attenuation problems–i.  _Nucleonics_ , 6(5):27, passim, 1950.
  * Kahn & Marshall (1953) Kahn, H. and Marshall, A. W.  Methods of reducing sample size in monte carlo computations.  _Journal of the Operations Research Society of America_ , 1(5):263–278, 1953.
  * Katharopoulos & Fleuret (2017) Katharopoulos, A. and Fleuret, F.  Biased importance sampling for deep neural network training.  _ArXiv_ , abs/1706.00043, 2017.  URL https://api.semanticscholar.org/CorpusID:38367260.
  * Katharopoulos & Fleuret (2018) Katharopoulos, A. and Fleuret, F.  Not all samples are created equal: Deep learning with importance sampling.  In Dy, J. and Krause, A. (eds.), _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 2525–2534. PMLR, 10–15 Jul 2018.  URL https://proceedings.mlr.press/v80/katharopoulos18a.html.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kondapaneni et al. (2019) Kondapaneni, I., Vévoda, P., Grittmann, P., Skřivan, T., Slusallek, P., and Křivánek, J.  Optimal multiple importance sampling.  _ACM Transactions on Graphics (TOG)_ , 38(4):37, 2019.
  * Krizhevsky et al. (2009) Krizhevsky, A., Hinton, G., et al.  Learning multiple layers of features from tiny images.  Technical report, Toronto, ON, Canada, 2009.
  * Langley (2000) Langley, P.  Crafting papers on machine learning.  In Langley, P. (ed.), _Proceedings of the 17th International Conference on Machine Learning (ICML 2000)_ , pp. 1207–1216, Stanford, CA, 2000. Morgan Kaufmann.
  * Loshchilov & Hutter (2015) Loshchilov, I. and Hutter, F.  Online batch selection for faster training of neural networks.  _arXiv preprint arXiv:1511.06343_ , 2015.
  * Needell et al. (2014) Needell, D., Ward, R., and Srebro, N.  Stochastic gradient descent, weighted sampling, and the randomized kaczmarz algorithm.  In Ghahramani, Z., Welling, M., Cortes, C., Lawrence, N., and Weinberger, K. (eds.), _Advances in Neural Information Processing Systems_ , volume 27. Curran Associates, Inc., 2014.  URL https://proceedings.neurips.cc/paper_files/paper/2014/file/f29c21d4897f78948b91f03172341b7b-Paper.pdf.
  * Owen & Zhou (2000) Owen, A. and Zhou, Y.  Safe and effective importance sampling.  _Journal of the American Statistical Association_ , 95(449):135–143, 2000.  doi: 10.1080/01621459.2000.10473909.  URL https://www.tandfonline.com/doi/abs/10.1080/01621459.2000.10473909.
  * Qi et al. (2017) Qi, C. R., Su, H., Mo, K., and Guibas, L. J.  Pointnet: Deep learning on point sets for 3d classification and segmentation.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 652–660, 2017.
  * Ren et al. (2019) Ren, H., Zhao, S., and Ermon, S.  Adaptive antithetic sampling for variance reduction.  In _International Conference on Machine Learning_ , pp. 5420–5428. PMLR, 2019.
  * Rumelhart et al. (1986) Rumelhart, D. E., Hinton, G. E., and Williams, R. J.  Learning representations by back-propagating errors.  _nature_ , 323(6088):533–536, 1986.
  * Santiago et al. (2021) Santiago, C., Barata, C., Sasdelli, M., Carneiro, G., and Nascimento, J. C.  Low: Training deep neural networks by learning optimal sample weights.  _Pattern Recognition_ , 110:107585, 2021.
  * Schaul et al. (2015) Schaul, T., Quan, J., Antonoglou, I., and Silver, D.  Prioritized experience replay.  _arXiv preprint arXiv:1511.05952_ , 2015.
  * Sitzmann et al. (2020) Sitzmann, V., Martel, J., Bergman, A., Lindell, D., and Wetzstein, G.  Implicit neural representations with periodic activation functions.  _Advances in neural information processing systems_ , 33:7462–7473, 2020.
  * Veach (1997) Veach, E.  _Robust Monte Carlo methods for light transport simulation_ , volume 1610.  Stanford University PhD thesis, 1997.
  * Wang et al. (2017) Wang, L., Yang, Y., Min, R., and Chakradhar, S.  Accelerating deep neural network training with inconsistent stochastic gradient descent.  _Neural Networks_ , 93:219–229, 2017.
  * Wu et al. (2015) Wu, Z., Song, S., Khosla, A., Yu, F., Zhang, L., Tang, X., and Xiao, J.  3d shapenets: A deep representation for volumetric shapes.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 1912–1920, 2015.
  * Zhang et al. (2019) Zhang, C., Öztireli, C., Mandt, S., and Salvi, G.  Active mini-batch sampling using repulsive point processes.  In _Proceedings of the AAAI conference on Artificial Intelligence_ , volume 33, pp. 5741–5748, 2019.
  * Zhang et al. (2023) Zhang, M., Dong, C., Fu, J., Zhou, T., Liang, J., Liu, J., Liu, B., Momma, M., Wang, B., Gao, Y., et al.  Adaselection: Accelerating deep learning training through data subsampling.  _arXiv preprint arXiv:2306.10728_ , 2023.
  * Zhao & Zhang (2015) Zhao, P. and Zhang, T.  Stochastic optimization with importance sampling for regularized loss minimization.  In Bach, F. and Blei, D. (eds.), _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pp. 1–9, Lille, France, 07–09 Jul 2015. PMLR.  URL https://proceedings.mlr.press/v37/zhaoa15.html.

##  Appendix A Optimization details

### Classification.

The classification tasks include image classification (MNIST (Deng, 2012),
CIFAR-10/100 (Krizhevsky et al., 2009) and point-cloud (ModelNet40 (Wu et al.,
2015)) classification.

The MNIST database contains 60,000 training images and 10,000 testing images.
We train a 3-layer fully-connected network (MLP) for MNIST over 50 epochs with
an Adam optimizer (Kingma & Ba, 2014). CIFAR-10, introduced by Krizhevsky et
al. (2009), is a dataset that consists of 60,000 color images of size 32x32.
These images belong to 10 different object classes, each class having 6,000
images. On the other hand, CIFAR-100 Krizhevsky et al. (2009) contains 100
classes with 600 images each. For each class, there are 500 training images
and 100 testing images. In our experiments, we train the ResNet-18 network (He
et al., 2016) on both datasets. We apply random horizontal flip and random
crops to augment the data during training. ModelNet40 contains 9,843 point
clouds for training and 2,468 for testing. Each point cloud has 1,024 points.
We train a PointNet (Qi et al., 2017) with 3 shared MLP layers and 2 fully-
connected layers for 300 epochs on point-cloud classification. We use the Adam
optimizer Kingma & Ba (2014), with batch size 64, weight decay 0.001, initial
learning rate 0.00002 divided by 2 after 100, 200 epochs.

### Regression.

Polynomial regression consists of optimizing the coefficients of a 1D
polynomial of a given order to fit randomly drawn data from a reference
polynomial of the same order. The reference data are generated on the interval
$[-2;2]$. Optimization is done using an Adam optimizer (Kingma & Ba, 2014)
with a mini-batch size of 32 elements.

The image regression task consists in learning the mapping between a 2D
coordinate input (pixel coordinate) and the 3-color output of the image for
this pixel. We use a network with 5 fully-connected layers associated with
positional encodings using SIREN activations (Sitzmann et al., 2020). The
training is done over 500 epoch using an Adam (Kingma & Ba, 2014) optimizer
and each mini-batch is composed of $256$ pixels for a $512^{2}$ reference
image.

##  Appendix B Multiple importance sampling in brief

### Importance sampling.

An Importance sampling Monte Carlo estimator $\langle F\rangle_{\mathrm{IS}}$
of a function $f$ is define as :

$\langle
F\rangle_{\mathrm{IS}}=\sum_{i=1}^{n}\frac{f(x_{i})}{np(x_{i})},\qquad
x_{i}\propto p(x).$ (9)

With $x_{i}$ the $i^{th}$ data sample drawn following the probability
distribution function $p(x)$.

The effectiveness of this estimator depends on the relation between the
functions $f(x)$ and $p(x)$. The variance of such estimator is :

$\mathrm{Var}{[\langle
F\rangle_{\mathrm{IS}}]}=\frac{1}{n}\mathrm{Var}{[\nicefrac{{f}}{{p}}]}.$ (10)

Reducing variance in the estimation depends on the proportionality between the
function $f$ and the probability density $p$.

When dealing with multivariate functions, finding a probability density
proportional to every parameters is often impractical. A trade-off is required
to obtain a single probability distribution maximizing the proportionality
with all the parameters of the function simultaneously. Several studies, such
as (Zhao & Zhang, 2015; Needell et al., 2014; Wang et al., 2017; Alain et al.,
2015), have shown that the optimal choice of sampling strategy is the $L_{2}$
norm of the function $f$.

### Multiple importance sampling.

Multiple Importance Sampling (MIS) is a technique that combines multiple
sampling strategies with associated weightings, unlike importance sampling
which relies on a single strategy. This approach allows for a more versatile
gradient estimation. The MIS Monte Carlo estimator, denoted as $\langle
F\rangle_{\mathrm{MIS}}$, is calculated by summing over all samples drawn
independently for each strategy, and then using a weighted estimator. The
equation for $\langle F\rangle_{\mathrm{MIS}}$ is given by:

$\langle
F\rangle_{\mathrm{MIS}}=\sum_{j=1}^{J}\sum_{i=1}^{n_{j}}w_{j}(x_{ij})\frac{f(x_{ij})}{n_{j}p_{j}(x_{ij})}$
(11)

Here, $x_{ij}$ represents the $i^{th}$ sample from the $j^{th}$ technique,
$w_{j}(x)$ is a weighting function such that $f(x)\neq
0\Rightarrow\sum^{J}_{j=1}w_{j}(x)=1$, and $p_{j}(x)=0\Rightarrow w_{j}(x)=0$.
$J$ is the number of sampling techniques, and $n_{j}$ is the number of samples
generated by the $j^{th}$ technique. The variance of a Monte Carlo estimator
using MIS, denoted as $\mathrm{Var}[\langle F\rangle_{\mathrm{MIS}}]$, can be
expressed as:

$\mathrm{Var}[\langle
F\rangle_{\mathrm{MIS}}]=\sum_{j=1}^{J}\int_{D}\frac{w_{j}(x)^{2}f(x)^{2}}{n_{j}p_{j}(x)}dx-\sum_{j=1}^{J}\frac{1}{n_{j}}\langle
w_{j},f\rangle^{2}$ (12)

The balance heuristic (Veach, 1997) is the most commonly used MIS heuristic.
It sets the weight of the samples from each technique according to the
following equation:

$w_{i}(x)=\frac{n_{i}p_{i}(x)}{\sum_{k=1}^{J}n_{k}p_{k}(x_{k})}$ (13)

This weighting strategy effectively mitigates the impact of events with low
probability when samples are drawn from a low-probability distribution. It
prevents a large increase in the contribution of such events in the Monte
Carlo estimator (11) where the function value would be divided by a very low
value. The balance heuristic compensates for this and avoids extreme cases.
Overall, this weighting strategy increases the robustness of the importance
sampling estimator, but it is limited by its independence from the function
value.

### Optimal weighting.

Following the discussion in Section 5, it can also be deduced from Equations
11 and 6 that $\langle F\rangle_{\mathrm{OMIS}}=\sum_{j=1}^{J}\alpha_{j}$.
Given a set of probability distribution functions $p_{1}$, …, $p_{J}$, we can
formulate the optimal MIS solver as Algorithm 3. $\boldsymbol{W}_{ij}$
represents the vector containing the balance weight (13) w.r.t. the J sampling
techniques and the normalization factor
$S(x_{ij})=\sum_{k=1}^{J}n_{k}p_{k}(x_{ij})$.

Algorithm 3 Optimal multiple importance sampling solver.

1:$\langle\boldsymbol{A}\rangle\leftarrow 0^{J\times
J},\langle\boldsymbol{b}\rangle\leftarrow 0^{J}$

2:for $t\leftarrow 1$ to $T$ do

3: for $j\leftarrow 1$ to $J$ do

4: $\\{x_{ij}\\}_{i=1}^{n_{j}}\leftarrow$ draw $n_{j}$ samples from technique
$p_{j}$

5:
$\langle\boldsymbol{A}\rangle\leftarrow\langle\boldsymbol{A}\rangle+\sum_{j=1}^{J}\sum_{i=1}^{n_{j}}\boldsymbol{W}_{ij}\boldsymbol{W}_{ij}^{T}$

6:
$\langle\boldsymbol{b}\rangle\leftarrow\langle\boldsymbol{b}\rangle+\sum_{j=1}^{J}\sum_{i=1}^{n_{j}}f(x_{ij})\boldsymbol{W}_{ij}/S(x_{ij})$

7:$\langle\boldsymbol{\alpha}\rangle\leftarrow\text{solve linear system
}\langle\boldsymbol{A}\rangle\langle\boldsymbol{\alpha}\rangle=\langle\boldsymbol{b}\rangle$

8:return $\sum_{j=1}^{N}\langle\boldsymbol{\alpha_{j}}\rangle$

The algorithm proceeds through three key stages. The first stage involves
initializing the linear system defined in Equation 7 (line 1). The second
stage iteratively updates the system for each drawn data sample (lines 5-6).
Upon completion of this process, the matrix $\boldsymbol{A}$ and vector
$\boldsymbol{b}$ provide Monte Carlo approximations of the quantities
specified in Equation 7. The third and final stage involve solving the linear
system to obtain the vector $\boldsymbol{\alpha}$ (line 7). The estimated
value of $\langle F\rangle_{\mathrm{MIS}}^{o}$ is then returned.

It can be noted that the linear system size scales with the number of sampling
techniques. More importantly each sampling technique needs to be sampled in
order create a linear system possible to solve. The number a sample of each
technique does not have to be the same but requires to be fixed at the start
of the algorithm. Also the presented algorithm works for a scalar value
function. In the case of multivariate function, multiple contribution vector
$\boldsymbol{b}$ need to be constructed (one per parameter) and the linear
system needs to be solved for each.

##  Appendix C Algorithm details

This section presents the two initialization subroutine for Algorithm 1 and
Algorithm 2. The role of the methods is to run a first epoch in a classical
SGD loop in order to process every data once. For each data the importance
metric is reported into the memory $q$ and returned with the current model
parameters. This approach avoids computing the importance for all data without
benefiting from the required forward step computed.

Algorithm 4 SGD-based initialization of persistent per-data importance $q$ in
Algorithm 1.

1:function Initialize($x$,$y$,$\Omega$,$\theta$,$B$)

2: for $t\leftarrow 1$ to $|\Omega|/B$ do

3: $x,y\leftarrow\\{x_{i},y_{i}\\}_{i=(t-1)\cdot{B}+1}^{t\cdot{B}+1}$

4: $l(x)\leftarrow{\mathcal{L}}(m(x,\theta),y)$

5: $\nabla l(x)\leftarrow$ Backpropagate$(l(x))$

6: $\langle\nabla L_{\theta}\rangle(x)\leftarrow$ $\nabla l(x)/B$ $\leftarrow$
Equation 3

7: $\theta\leftarrow\theta-\eta\,\langle\nabla L_{\theta}\rangle(x)$
$\leftarrow$ Equation 2

8: $q(x)\leftarrow\left\|\frac{\partial\mathcal{L}(x)}{\partial
m(x,\theta)}\right\|$

9: return $q$,$\theta$

Algorithm 5 Subroutine for initialization in Algorithm 2

1:function InitializeMIS($x$,$y$,$\Omega$,$\theta$,$B$)

2: Initialize $\boldsymbol{q}$ in a classical SGD loop

3: for $t\leftarrow 1$ to $|\Omega|/B$ do

4: $x,y\leftarrow\\{x_{i},y_{i}\\}_{i=(t-1)\cdot{B}+1}^{t\cdot{B}+1}$

5: See all samples in the first epoch

6: $l(x)\leftarrow{\mathcal{L}}(m(x,\theta),y)$

7: $\nabla l(x)\leftarrow$ Backpropagate$(l(x))$

8: $\langle\nabla L_{\theta}\rangle(x)\leftarrow$ $\nabla l(x)/B$ $\leftarrow$
Equation 3

9: $\theta\leftarrow\theta-\eta\,\langle\nabla L_{\theta}\rangle(x)$
$\leftarrow$ Equation 2

10: $\boldsymbol{q}(x)\leftarrow\frac{\partial\mathcal{L}(x)}{\partial
m(x,\theta)}$

11: return $\boldsymbol{q}$,$\theta$

Algorithm 6 Subroutine for cross entropy loss importance metric

1:$x_{i}=$ data sample, $y_{i}=$ class index of $x_{i}$

2:function Importance($x_{i}$,$y_{i}$)

3: $s\leftarrow\exp(m(x_{i},\theta))/\sum_{k=1}^{C}\exp(m(x_{i},\theta)_{k})$
$\leftarrow$ Eq.14

4: $q\leftarrow\sum_{j=1}^{C}s_{j}-\mathbf{1}_{j=y_{i}}$ $\leftarrow$ Eq.16

5: return $q$

##  Appendix D Cross-entropy loss gradient

Machine learning frameworks take data $x$ as input, perform matrix
multiplication with weights and biases added. The output layer is then fed to
the softmax function to obtain values ${s}$ that are fed to the loss function.
$y$ represents the target values. We focus on the categorical cross-entropy
loss function for the classification problem (with $C$ categories) given by:

${\mathcal{L}}_{\text{cross-
ent}}=-\sum_{i}y_{i}\log{s}_{i},\;\;\;{s}_{i}=\frac{\exp(m(x_{i},\theta)_{l})}{\sum_{l}^{C}\exp(m(x_{i},\theta)_{l})}.\\!$
(14)

For backpropagation, we need to calculate the derivative of the $\log{s}$ term
w.r.t. the weighted input $z$ of the output layer. We can easily derive the
derivative of the loss from first principles as shown below:

$\begin{split}\frac{\partial{\mathcal{L}}_{\text{cross-ent}}}{\partial
m(x_{i},\theta)_{j}}&=-\frac{\partial}{\partial
m(x_{i},\theta)_{j}}\left(\sum_{i}^{C}y_{i}\log{s}_{i}\right)\\\
&=-\sum_{i}^{C}y_{i}\frac{\partial}{\partial
m(x_{i},\theta)_{j}}\log{s}_{i}\\\
&=-\sum_{i}^{C}\frac{y_{i}}{{s}_{i}}\frac{\partial{s}_{i}}{\partial
m(x_{i},\theta)_{j}}\\\
&=-\sum_{i}^{C}\frac{y_{i}}{{s}_{i}}{s}_{i}\cdot(\mathbf{1}\\{i==j\\}-{s}_{j})\\\
&=\sum_{i}^{C}{y_{i}}\cdot{s}_{j}-\sum_{i}^{C}y_{i}\cdot(\mathbf{1}\\{i==j\\})\\\
&={s}_{j}\sum_{i}^{C}{y_{i}}-y_{j}={s}_{j}-y_{j}\end{split}$ (15)

The partial derivative of the cross-entropy loss function w.r.t. output layer
parameters has the form:

$\displaystyle\frac{\partial{\mathcal{L}}_{\text{cross-ent}}}{\partial
m(x_{i},\theta)_{j}}$ $\displaystyle={s}_{j}-y_{j}$ (16)

For classification tasks, we directly use this analytic form of the derivative
and compute it’s norm as weights for adaptive and importance sampling.

##  Appendix E Importance momentum

Updating the persistent per-sample importance $q$ directly sometime leads to a
sudden decrease of accuracy during training. To make the training process more
stable, we update $q$ by linearly interpolating the importance at the previous
and current steps:

$q(x)=\gamma\cdot q_{prev}(x)+(1-\gamma)\cdot q(x)$ (17)

where $\gamma$ is a constant for all data samples. In practice, we use
$\gamma\in\\{0.0,0.1,0.2,0.3\\}$ as it gives the best trade-off between
importance update and stability. This can be seen as a momentum evolution of
the per-sample importance to avoid high variation. Utilizing an exponential
moving average to update the importance metric prevents the incorporation of
outlier values. This is particularly beneficial in noisy setups, like
situations with a high number of class or a low total number of data.

##  Appendix F Additional results

This section provides additional results, including an ablation study as shown
in Figure 9 for linear-system momentum used in Algorithm 2 and results of our
adaptive sampling method. Figures 9 and 6 demonstrate that classical SGD, DLIS
and Our IS work similarly with and without momentum. Our OMIS outperforms
other methods in both cases.

Figures 8, 10 and 11 show that our adaptive sampling variant (our AS) can
achieve better results than our IS or our OMIS in practice. Our AS is a
heuristic and we leave its theoretical formulation as future work.

Figure 8: We compare loss for the MNIST dataset between the resampling
algorithm by Katharopoulos & Fleuret (2018) (DLIS) and our algorithm. At equal
epochs, DLIS works better than both classical and resampling SGD. However, at
equal time, the resampling cost is too high, making DLIS even slower than
standard SGD.  Figure 9:  Ablation study on point-cloud classification using
linear-system momentum as described in Algorithm 2 for baselines represented
as dashed lines. Our OMIS still outperforms other baselines at equal epochs,
similar to the results shown in Figure 6.  Figure 10: Comparisons on CIFAR-10
using Vision Transformer (ViT) (Dosovitskiy et al., 2020). The results show
our importance sampling scheme (Our IS) and the adaptive sampling variant (Our
AS) can improve over classical SGD, LOW (Santiago et al., 2021) and DLIS
(Katharopoulos & Fleuret, 2018) on modern transformer architecture.  Figure
11:  On CIFAR-100 classification dataset, instead of comparing the DLIS
resampling algorithm, we use DLIS importance in our Algorithm 1. We display
zoom-in of the end of the curves to highlight the differences. At equal epochs
(left), our methods (Our IS & Our AS) show improvements compared to LOW
(Santiago et al., 2021) and DLIS weights. At equal time (right), LOW and the
DLIS weights takes longer to converge. Overall our approach shows faster
convergence with lower importance computation.