# Theoretical model for en face optical coherence tomography imaging and its
application to volumetric differential contrast imaging

Kiriko Tomita  1 Shuichi Makita  1 Naoki Fukutake  2 Rion Morishita  1 Ibrahim
Abd El-Sadek  1,3 Pradipta Mukherjee  1 Antonia Lichtenegger  1,4 Junya
Tamaoki  5 Lixuan Bian  5 Makoto Kobayashi  5 Tomoko Mori  6 Satoshi Matsusaka
6 and Yoshiaki Yasuno1  1Computational Optics Group, University of Tsukuba,
Japan
2Nikon Corporation, Japan
3Department of Physics, Faculty of Science, Damietta University, Damietta,
Egypt
4Center for Medical Physics and Biomedical Engineering, Medical University of
Vienna, Austria
5Department of Molecular and Developmental Biology, Institute of Medicine,
University of Tsukuba, Japan
6Clinical Research and Regional Innovation, Faculty of Medicine, University of
Tsukuba, Japan
<EMAIL_ADDRESS>
###### Abstract

A new formulation of lateral imaging process of point-scanning optical
coherence tomography (OCT) and a new differential contrast method designed by
using this formulation are presented. The formulation is based on a
mathematical sample model called the dispersed scatterer model (DSM), in which
the sample is represented as a material with a spatially slowly varying
refractive index and randomly distributed scatterers embedded in the material.
It is shown that the formulation represents a meaningful OCT image and speckle
as two independent mathematical quantities. The new differential contrast
method is based on complex signal processing of OCT images, and the physical
and numerical imaging processes of this method are jointly formulated using
the same theoretical strategy as in the case of OCT. The formula shows that
the method provides a spatially differential image of the sample structure.
This differential imaging method is validated by measuring in vivo and in
vitro samples.

††journal: boe††articletype: Research Article

https://optics.bk.tsukuba.ac.jp/COG/

##  1 Introduction

Several theoretical models of optical coherence tomography (OCT) [1] have been
studied to understand its image formation. Because OCT is used in the
biomedical field such as ophthalmology [2] and dermatology [3], theoretical
modeling is important for interpreting OCT images and clarifying what is
visualized in the sample tissue. Theoretical modeling is also useful for
designing optimal OCT systems and developing image processing algorithms.

Several analytic and numerical models of OCT imaging have been developed.
After an early trial of OCT modeling [4], the extended Huygens-Fresnel
formulation was introduced [5, 6]. In this model, propagation of the OCT probe
beam in a tissue, i.e., an inhomogeneous medium, is modeled by extended
Huygens-Fresnel formulation. Numerical simulations based on this model can
recapitulate the forward multiple scatting in the tissue. Monte-Carlo methods
were also widely used to numerically emulate the OCT imaging [7, 8]. Although
Monte-Carlo methods are nondeterministic and computationally expensive in
general, they can be applied to wide varieties of imaging models, conditions,
and samples. The interaction between the probe beam and the sample can be also
accurately computed by using Maxwell’s equation [9, 10]. An extension of this
method, the pseudo-spectral time domain method, enables the realistic
simulation of three-dimensional OCT imaging [11]. Recently, the speckle was
mathematically described as the interference of light backscattered by
multiple scatterers in a resolution volume [12].

Aforementioned approaches aimed to recapture accurate and realistic OCT
signals. By contrast, we aim to establish a simplified theoretical framework
that enables the intuitive understanding of the point-scanning OCT imaging
process and the interpretation of an en face OCT image. In this framework, we
model a sample as a spatially slowly varying refractive index in which
scatterers are randomly dispersed (Section 2). The density of the scatterers
is also spatially slowly varying. This representation of the sample is denoted
as a “dispersed scatterer model” (DSM). We formulate the OCT imaging process
using the DSM (Section 3). This formulation identifies the relationship
between the sample structure, that is, the slowly varying refractive index
distribution and the scatterer density distribution, and the en face OCT
signal. In this formulation, a meaningful OCT signal and speckle appear as
different terms; hence, it provides a comprehensive understanding of the en
face OCT image formation process.

By actively using this theoretical framework, we design and experimentally
demonstrate an OCT-based volumetric differential contrast (VDC) imaging method
(Section 4). We implement this method as the numerical signal processing flow
of complex en face OCT images and provide differential contrast images similar
to those of differential interference contrast microscopy (DIC) [13]. Similar
to DIC, VDC highlights the micro-structure of a sample, but unlike DIC, VDC
can be volumetric.

The DSM-based framework can also explain phenomena that occur in OCT-based
imaging methods but have not been theoretically well explained. An example is
artifacts found in A-scan-wise Doppler OCT images caused by the structure of
superior tissue. Another example relates to the specific imaging condition
required for point-scanning computational multi-directional OCT [14]. In
experiments, it was found that defocus is mandatory to achieve this method
[15], but this necessity has not been explained in theory. We analyze these
issues using the DSM-based framework in Section 5.

##  2 Dispersed scatterer model (DSM)

In our theoretical framework, we model a sample to be measured as a spatially
slowly varying refractive index distribution in which scatterers are spatially
randomly dispersed. We consider the density of the scatterers to be also
spatially slowly varying. We refer to this theoretical model of the sample as
a DSM.

Figure 1 schematically illustrates the DSM, where $(x,y,z)$ is the position in
the sample as $x$ and $y$ are the lateral positions, and $z$ is the depth
position. The spatially slowly varying refractive index is represented by
$n(x,y,z)$, in which small scatterers are dispersed (i.e., embedded), where
the scatterers are assumed to be significantly smaller than the spatial
variation of $n(x,y,z)$. The position of the $i$-th scatterer is
$(x_{i},y_{i},z_{i})$ and all the scatterers have an identical refractive
index of $n_{s}$. The probe beam is incident from above of the sample, and the
green region in the figure represents the coherence gate. In our formulations
for OCT (Section 3) and VDC (Section 4), we consider en face images in which
scattering only within this coherence gate contributes to imaging.

Figure 1:  The schematic model of the sample, i.e., the dispersed scatterer
model (DSM). The sample is represented as a spatially slowly varying
refractive index $n(x,y,z)$ in which infinitely small scatterers (black dots)
are randomly dispersed. All the scatterers are assumed to have the same
refractive index of $n_{s}$. The green region represents the coherence gate of
OCT and only the scatterers within the coherence gate contribute to OCT
imaging. The probe beam is incident on the sample from the top and only lights
backscattered toward the incident direction contribute to imaging.

##  3 Formulation of OCT imaging

In this section, we formulate the image formation process of standard scanning
OCT imaging using DSM.

###  3.1 Complex OCT signal

We assume that a Gaussian beam is incident on the sample. Considering that $N$
scatterers are involved in imaging, i.e., $N$ scatterers exist in the
coherence gate, a complex en face OCT signal becomes

$s(\bm{x^{\prime}};z^{\prime}_{0})=\sum_{i=1}^{N}B({\bm{x}_{i}},z_{0})e^{i\phi_{i}}{P_{c}\left(\bm{x^{\prime}}-{\bm{x}_{i}};z^{\prime}_{0}\right)},$
(1)

where $\bm{x}=(x,y)$ and $\bm{x^{\prime}}=(x^{\prime},y^{\prime})$ are the
two-dimensional lateral positions on a sample and image, respectively. Here
the contribution from each scatterer is that of single scattering. Namely, the
first-order Born approximation has been employed. Throughout this paper, we
use $\bm{x}$ and $z$ to represent the position in the object (sample) space,
and $\bm{x^{\prime}}$ and $z^{\prime}$ to represent the position in the image
space. $z_{0}^{\prime}$ is the depth position in the image corresponding to
the center depth of coherence gate $z_{0}$ in the sample. The magnification
between the sample and the image is assumed to be 1.
${\bm{x}_{i}}=(x_{i},y_{i})$ is the lateral position of the $i-$th scatterer.
$B(\bm{x},z_{0})$ is a quantity defined by the slowly varying refractive index
distribution of the medium, $n(x,y,z)$, and the refractive index of the
embedded scatterers, $n_{s}$. And hence, $B$ is also spatially slowly varying
as is $n$. It is noteworthy that the overall profile of $B(\bm{x},z_{0})$ is
defined solely by the profile of the spatially slowly varying refractive index
distribution $n(x,y,z)$. The spatial inhomogeneity of the dispersed scatterers
is accounted for later when we derive the continuous form of the imaging
formula in Section 3.3. For the case of sufficiently small scatterers,
$B(\bm{x},z_{0})$ can be considered as the scattering potential
$B(\bm{x},z_{0})\simeq{\left\\{{2\pi
n(\bm{x},z_{0})/\lambda_{0}}\right\\}}\left[\left\\{n_{s}/n(\bm{x},z_{0})\right\\}^{2}-1\right]$,
where $\lambda_{0}$ is the central wavelength of the probe beam [16].

${P_{c}\left(\bm{x};z^{\prime}_{0}\right)}$ is a complex point spread function
(PSF) at depth $z^{\prime}_{0}$. $\phi_{i}$ is the phase offset caused by the
depth position of the $i$-th scatterer with respect to the central depth of
the coherence gate $z_{0}$:

$\phi_{i}=k_{0}\,n\,2(z_{i}-z_{0}),$ (2)

where $k_{0}=2\pi/\lambda_{0}$ is the wavenumber as $\lambda_{0}$ is the
center wavelength of the probe beam and $n$ is the refractive index of the
sample around the scatterers. Note that we assume that $B$ is homogeneous
along the depth within the coherence gate; hence, the variation of the
scatterer’s depth position does not affect $B$ but only appears as phase
offset.

The complex PSF ${P_{c}\left(\bm{x};z^{\prime}_{0}\right)}$ is further broken
down as

${P_{c}\left(\bm{x};z^{\prime}_{0}\right)}\equiv{{P_{a}}\left(\bm{x};z^{\prime}_{0}\right)}\,e^{i\varphi\left(\bm{x};z_{0}\right)},$
(3)

where $P_{a}\equiv{\left|{P_{c}}\right|}$ is the amplitude of the PSF and is
typically Gaussian in the lateral direction, and $\varphi$ is the phase of the
complex PSF. We discuss a more practical shape for the complex PSF in Section
3.4.

By substituting Eq. (3) into Eq. (1), the complex OCT signal
${s\left(\bm{x^{\prime}};z^{\prime}_{0}\right)}$ becomes

${s\left(\bm{x^{\prime}};z^{\prime}_{0}\right)}=\sum_{i=1}^{N}B({\bm{x}_{i}},z_{0})e^{i\phi_{i}}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{i}};z^{\prime}_{0}\right)}\,e^{i\varphi\left(\bm{x^{\prime}}-{\bm{x}_{i}};z_{0}\right)}.$
(4)

As shown in Eq. (4), the contribution from each scatterer is weighted by the
amplitude of the complex PSF, $P_{a}$, and its phase is modulated by both the
depth position of the scatterer and spatial phase pattern of the PSF.

###  3.2 Intensity-OCT signal in discrete form

The intensity of OCT is the absolute square of the complex OCT signal, that
is, the product of the complex OCT signal and its conjugate. From Eq. (4), the
en face OCT intensity $I\left(\bm{x^{\prime}};z_{0}\right)$ becomes

$\displaystyle I\left(\bm{x^{\prime}};z^{\prime}_{0}\right)$
$\displaystyle={s\left(\bm{x^{\prime}};z^{\prime}_{0}\right)}{s^{*}\left(\bm{x^{\prime}};z^{\prime}_{0}\right)}$
$\displaystyle=\left[\sum_{i=1}^{N}B({\bm{x}_{i}},z_{0})e^{i\phi_{i}}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{i}};z^{\prime}_{0}\right)}\,e^{i\varphi\left(\bm{x^{\prime}}-{\bm{x}_{i}}\right)}\right]$
$\displaystyle\quad\times\left[\sum_{j=1}^{N}B({\bm{x}_{j}},z_{0})e^{-i\phi_{j}}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{j}};z^{\prime}_{0}\right)}\,e^{-i\varphi\left(\bm{x^{\prime}}-{\bm{x}_{j}}\right)}\right]$
$\displaystyle=\sum_{i=1}^{N}\sum_{j=1}^{N}B({\bm{x}_{i}})B({\bm{x}_{j}})e^{i(\phi_{i}-\phi_{j})}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{i}}\right)}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{j}}\right)}e^{i\left[\varphi\left(\bm{x^{\prime}}-{\bm{x}_{i}}\right)-\varphi\left(\bm{x^{\prime}}-{\bm{x}_{j}}\right)\right]},$
(5)

where ${\bm{x}_{j}}=(x_{j},y_{j})$ is the position of the $j$-th scatterer in
the sample. On the last line of the equation, we did not explicitly write
$z_{0}$ and $z^{\prime}_{0}$ for simplicity. We also omitted them in the later
part of the paper if it did not cause confusion.

Figure 2:  (a) Schematic of the formulation of OCT intensity. Each circle
represents each term of the equation. The same color indicates the same term
and the asterisk represents the complex conjugate. The red box consists of the
terms on the second line of Eq. (5) and the blue box corresponds to the terms
on the third line of Eq. (5). The orange box represents the terms on the
fourth line of Eq. (5). (b) Physical depiction of the numerical interactions
between two signals contributed by each pair of scatterers, which are
mathematically described by Eq. (6). The interactions are classified into two
groups: auto-interaction (AI) and cross-interaction (CI) groups. AI
corresponds to the numerical interference of the electric field scattered from
a scatterer and the conjugate field from the same scatterer; hence, it
corresponds to the backscattered intensity from the scatterer. CI corresponds
to the numerical interference of the electric field from a scatterer and the
conjugate field from another scatterer; hence, is a phasor. AI and CI are also
referred to as auto-interference and cross-interference, respectively.

Eq. (5) is schematically depicted in Fig. 2(a). Each circle represents each
term of the equation, that is, the contribution from each scatterer in the
sample. The same color indicates the same term, that is, the contribution from
the same scatterer, and the asterisk represents the complex conjugate. The red
box consists of the terms in ${s\left(\bm{x^{\prime}},z^{\prime}_{0}\right)}$,
that is, the second line of Eq. (5) and the blue box represents
${s^{*}\left(\bm{x^{\prime}},z^{\prime}_{0}\right)}$, which corresponds to the
third line of Eq. (5). The orange box represents interactions among the terms
in ${s\left(\bm{x^{\prime}},z^{\prime}_{0}\right)}$ and
${s^{*}\left(\bm{x^{\prime}},z^{\prime}_{0}\right)}$, that is, the terms on
the fourth line of Eq. (5). We use this depiction in the interpretation of the
equation in the next paragraph.

The terms of Eq. (5) can be classified into two groups for $i=j$ and $i\neq j$
as

$I\left(\bm{x^{\prime}}\right)=\color[rgb]{0.2,0.2,0.6}\sum_{i=1}^{N}B^{2}({\bm{x}_{i}})P_{a}^{2}\left(\bm{x^{\prime}}-{\bm{x}_{i}}\right)\\\
+\color[rgb]{0.43,0.21,0.1}\underset{(i\neq
j)}{\sum_{i=1}^{N}\sum_{j=1}^{N}}B({\bm{x}_{i}})B({\bm{x}_{j}})\color[rgb]{0.8,0.0,0.0}e^{i(\phi_{i}-\phi_{j})}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{i}}\right)}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{j}}\right)}e^{i\left[\varphi\left(\bm{x^{\prime}}-{\bm{x}_{i}}\right)-\varphi\left(\bm{x^{\prime}}-{\bm{x}_{j}}\right)\right]}\color[rgb]{0,0,0}.$
(6)

This equation is the representation of the OCT intensity signal in the
discrete form of our theoretical framework. In this formulation, each term in
the first group (the first line of the equation, for $i=j$, written in blue)
is the numerical interaction between contributions from the same scatterer;
hence, we refer to this group as an auto-interaction (AI) group. By contrast,
each term in the second group (the second line, for $i\neq j$, written in
brown and red) is the numerical interaction between contributions from
different scatterers. We refer to this group as a cross-interaction (CI)
group. We classified CI into brown and red parts, and this classification is
used in later sections. The terms of the AI group correspond to the diagonal
terms in the orange box in Fig. 2(a), whereas those of the CI group are the
off-diagonal terms.

It is noteworthy that these interactions can be considered as the numerical
interferences of scattered electric fields from the scatterers as depicted in
Fig. 2(b). The AI term corresponds to interference between an electric field
from a scatterer and the complex conjugate of the same electric field. Hence,
each term in the AI group represents the squared intensity of scattered light
from each scatterer and the AI group is the incoherence summation of these
scattered light intensities. By contrast, each term in the CI group is the
numerical interference between the scattered complex field from a scatterer
and the conjugate complex field from another scatterer, where two scatterers
are sufficiently close to be included in the complex PSF. Note that each term
in the CI group is a phasor with a practically random phase and each term in
the CI group has its complex conjugate pair, as can be schematically seen from
the diagonal symmetricity of the orange box in Fig. 2(a). Because of the
random phase of each term, the summation of a term and its complex conjugate
becomes a random real value. This suggests that the CI group causes a random
signal, i.e., speckle, as we discuss in detail in Section 3.5.

###  3.3 Intensity-OCT signal in continuous form

In Eq. (6), the OCT intensity is expressed using the contributions from
discrete scatterers; hence, the equation is also in discrete form. In this
subsection, we derive a continuous-form representation of the OCT intensity,
which is a convolution-based conventional image formation equation, and hence
convenient to interpret.

Using the Dirac delta function, the intensity-OCT signal in discrete form [Eq.
(6)] can be rewritten as

$I(\bm{x^{\prime}})=\color[rgb]{0.2,0.2,0.6}\iint\sum_{i=1}^{N}\delta{\left({{\bm{x}_{p}}-{\bm{x}_{i}}}\right)}B^{2}({\bm{x}_{p}})P_{a}^{2}(\bm{x^{\prime}}-{\bm{x}_{p}})d{\bm{x}_{p}}\\\
+\color[rgb]{0.43,0.21,0.1}\biggl{[}\iiiint\sum_{i=1}^{N}\sum_{j=1}^{N}\delta{\left({{\bm{x}_{p}}-{\bm{x}_{i}},{\bm{x}_{q}}-{\bm{x}_{j}}}\right)}B({\bm{x}_{p}})B({\bm{x}_{q}})\color[rgb]{0.8,0.0,0.0}e^{i(\phi_{pq})}\\\
\color[rgb]{0.8,0.0,0.0}\times{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{p}}\right)}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{q}}\right)}e^{i\left[\varphi\left(\bm{x^{\prime}}-{\bm{x}_{p}}\right)-\varphi\left(\bm{x^{\prime}}-{\bm{x}_{q}}\right)\right]}\color[rgb]{0.43,0.21,0.1}d{\bm{x}_{p}}d{\bm{x}_{q}}\\\
\color[rgb]{0.43,0.21,0.1}-\iint\sum_{i=1}^{N}\delta{\left({{\bm{x}_{p}}-{\bm{x}_{i}}}\right)}B^{2}({\bm{x}_{p}})P_{a}^{2}(\bm{x^{\prime}}-{\bm{x}_{p}})d{\bm{x}_{p}}\biggr{]}\color[rgb]{0,0,0},$
(7)

where ${\bm{x}_{p}}$ and ${\bm{x}_{q}}$ are variables in real physical space
in the sample. Each of the double integral is for $x$ and $y$ components of
the vectorized integration variable ${\bm{x}_{p}}$, while the quadruple
integral is for two vectorized integration variables ${\bm{x}_{p}}$ and
${\bm{x}_{q}}$. The intervals of all the integrals are from $-\infty$ to
$+\infty$. Hereafter, the integral intervals of all equations are from
$-\infty$ to $+\infty$; we omit them for simplicity. $\phi_{pq}$ is a phase
that is originated from $\phi_{i}-\phi_{j}$ in Eq. (6). Namely, $\phi_{pq}$ is
the phase difference between the contributions of $i$-th and $j$-th scatterer
and equals to $\phi_{i}-\phi_{j}$. Similar to $\phi_{i}-\phi_{j}$, $\phi_{pq}$
is unpredictable and random in practice, and $\phi_{pq}=-\phi_{qp}$. Note
that, in the second line of Eq. (7), we can consider that ${\bm{x}_{p}}$ and
${\bm{x}_{q}}$ respectively denote the positions of $i$-th and $j$-th
scatterers due to the Dirac delta function. The first line of Eq. (7) (blue)
corresponds to the AI group [the first line of Eq. (6)] and the second to
fourth lines of Eq. (7) (brown and red) correspond to the CI group [the second
line of Eq. (6)]. Note that the summation in CI of Eq. (6) is for $i\neq j$,
whereas we removed this condition in Eq. (7). Alternatively, we subtract the
fourth line, which is identical to the first line corresponding to the case of
$i=j$. This removal of condition was necessary to convert the discrete
equation [Eq. (6)] into the continuous equation [Eq. (7)] because the integral
cannot properly account for such a condition.

We introduce functions $D$ and $D^{\prime}$ as

$D{\left({{\bm{x}_{p}}}\right)}\equiv\sum_{i=1}^{N}\delta{\left({{\bm{x}_{p}}-{\bm{x}_{i}}}\right)},$
(8)
$D^{\prime}{\left({{\bm{x}_{p}},{\bm{x}_{q}}}\right)}\equiv\sum_{i=1}^{N}\sum_{j=1}^{N}\delta{\left({{\bm{x}_{p}}-{\bm{x}_{i}},{\bm{x}_{q}}-{\bm{x}_{j}}}\right)},$
(9)

which are two- and four-dimensional binary maps of the scatterers’ positions,
respectively. When the scatterer density is high enough to be considered as
continuously distributed, that is, $N$ is sufficiently large, $D$ and
$D^{\prime}$ can be considered as non-binary scatterer density maps in two-
and four-dimensional spaces, respectively. Specifically, $D$ is the scatterer
density map of the sample.

By substituting Eqs. (8) and (9) into Eq. (7), we obtain the continuous form
of the equation:

$I{\left({\bm{x^{\prime}}}\right)}=\color[rgb]{0.2,0.2,0.6}\iint
D({\bm{x}_{p}})B^{2}({\bm{x}_{p}})P_{a}^{2}{\left({\bm{x^{\prime}}-{\bm{x}_{p}}}\right)}d{\bm{x}_{p}}\\\
+\color[rgb]{0.43,0.21,0.1}\left[\iiiint
D^{\prime}({\bm{x}_{p}},{\bm{x}_{q}})B({\bm{x}_{p}})B({\bm{x}_{q}})\right.\\\
\color[rgb]{0.8,0.0,0.0}\left.e^{i(\phi_{pq})}\right.\color[rgb]{0.8,0.0,0.0}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{p}}\right)}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{q}}\right)}e^{i\left[\varphi\left(\bm{x^{\prime}}-{\bm{x}_{p}}\right)-\varphi\left(\bm{x^{\prime}}-{\bm{x}_{q}}\right)\right]}\color[rgb]{0.43,0.21,0.1}d{\bm{x}_{p}}d{\bm{x}_{q}}\\\
\color[rgb]{0.43,0.21,0.1}\left.-\iint
D({\bm{x}_{p}})B^{2}({\bm{x}_{p}})P_{a}^{2}{\left({\bm{x^{\prime}}-{\bm{x}_{p}}}\right)}d{\bm{x}_{p}}\right]\color[rgb]{0,0,0}.$
(10)

The first integral (i.e., the AI part) of Eq. (10) (blue) is “an imaging
equation of the intensity-OCT image” in continuous form. According to this
equation, the object to be imaged in the intensity-OCT image is the product of
the squared scattering potential $B^{2}$ and the scatterer density $D$, and
the PSF of intensity-OCT imaging is $P_{a}^{2}$. A more detailed
interpretation of this intensity-OCT imaging formula is provided in the
following subsections (Sections 3.4 and 3.5).

###  3.4 Intensity-OCT imaging with Gaussian PSF

The property of intensity-OCT imaging becomes clearer if we assume a complex
Gaussian PSF with defocus in the continuous OCT imaging formula [Eq. (10)].
The complex Gaussian PSF is obtained from a Gaussian probe beam, and the
Gaussian probe beam is typical for a single-mode-fiber-based point-scanning
OCT. According to the formulation by Ralston et al.[17], the Gaussian probe
beam in the paraxial approximation becomes

$G(\bm{x},z;z_{0})=\frac{w_{0}}{w(z;z_{0})}\exp\left[{-\frac{\bm{x}^{2}}{w^{2}(z;z_{0})}}\right]\exp\left[-i\left\\{nk_{0}(z-z_{0})+\frac{nk_{0}\bm{x}^{2}}{2R(z;z_{0})}-\psi(z;z_{0})\right\\}\right],$
(11)

where $z_{0}$ is the depth position of the focus, $w(z;z_{0})$ is the radius
of the probe beam at depth $z$, $w_{0}=w(z=z_{0};z_{0})$ is the beam radius at
the in-focus depth, that is, the beam waist radius, and $n$ is the refractive
index of the sample around the position of interest. The second term in the
second exponential represents the quadratic phase induced by defocus and
$R(z;z_{0})$ is the phase curvature. $\bm{x}^{2}$ represents
$\bm{x}\cdot\bm{x}$. Here after the squared vector represents the inner
product of a vector with itself. In the third term of the second exponential,
$\psi(z;z_{0})$ represents the Gouy phase.

The beam waist radius $w_{0}$, the radius of the probe beam $w$, and the phase
curvature $R$ are

$w_{0}=\frac{4f}{n\,\Phi\,k_{0}},$ (12)
$w(z;z_{0})=w_{0}\sqrt{1+{\left\\{{\frac{2(z-z_{0})}{n^{2}\,k_{0}\,w_{0}^{2}}}\right\\}}^{2}},$
(13)
$R(z;z_{0})={\left({z-z_{0}}\right)}{\left[{1+{\left\\{{\frac{n^{2}\,k_{0}\,w_{0}^{2}}{2(z-z_{0})}}\right\\}}^{2}}\right]},$
(14)

where $f$ is the focal length of the objective and $\Phi$ is the
$1/e^{2}$-diameter of the probe beam incident on the objective.

Because OCT is a reflective imaging modality, and the illumination and
collection paths share the same optics, the complex PSF
${P_{c}\left(\bm{x},z;z_{0}\right)}$ is obtained by multiplying
$G(\bm{x},z;z_{0})$ by itself as

${P_{c}\left(\bm{x},z;z_{0}\right)}=G(\bm{x},z;z_{0})\,G(\bm{x},z;z_{0})\\\
=\frac{w_{0}^{2}}{w^{2}(z;z_{0})}\exp\left[-\frac{2\bm{x}^{2}}{w^{2}(z;z_{0})}\right]\exp\left[-i\left\\{2nk_{0}(z-z_{0})+\frac{nk_{0}\bm{x}^{2}}{R(z;z_{0})}-2\psi(z;z_{0})\right\\}\right].$
(15)

On the first line, each $G(\bm{x},z;z_{0})$ is for probe beam illumination and
collection, respectively.

By comparing this complex PSF and that of Eq. (3), we find that

${{P_{a}}\left(\bm{x},z^{\prime};z^{\prime}_{0}\right)}=\frac{w_{0}^{2}}{w^{2}(z^{\prime};z^{\prime}_{0})}\exp\left[-\frac{2\bm{x}^{2}}{w^{2}(z^{\prime};z^{\prime}_{0})}\right]$
(16)

and

$\varphi(\bm{x},z;z_{0})=-2nk_{0}(z-z_{0})-\frac{nk_{0}}{R(z;z_{0})}\bm{x}^{2}+2\psi(z;z_{0}).$
(17)

By substituting this Gaussian complex PSF, or equally $P_{a}$ and $\varphi$,
into Eq. (6), we obtain the discrete-form intensity-OCT imaging equation:

$I(\bm{x^{\prime}})=\color[rgb]{0.2,0.2,0.6}\frac{w_{0}^{4}}{w^{4}}\sum_{i=1}^{N}B^{2}({\bm{x}_{i}})\exp{\left[{-\frac{4}{w^{2}}\left(\bm{x^{\prime}}-{\bm{x}_{i}}\right)^{2}}\right]}\\\
+\color[rgb]{0.43,0.21,0.1}\frac{w_{0}^{4}}{w^{4}}\underset{(i\neq
j)}{\sum_{i=1}^{N}\sum_{j=1}^{N}}B({\bm{x}_{i}})B({\bm{x}_{j}})\color[rgb]{0.8,0.0,0.0}\exp{\left[{-\frac{({\bm{x}_{i}}-{\bm{x}_{j}})^{2}}{w^{2}}}\right]}\exp
i{\left[{\phi_{i}-\phi_{j}}\right]}\\\
\color[rgb]{0.8,0.0,0.0}\times\exp{\left[{-\frac{4}{w^{2}}{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{i}}+{\bm{x}_{j}}}{2}}\right)}^{2}}\right]}\exp
i{\left[{\frac{2nk_{0}}{R}{\left({{\bm{x}_{i}}-{\bm{x}_{j}}}\right)}\cdot{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{i}}+{\bm{x}_{j}}}{2}}\right)}}\right]}\color[rgb]{0,0,0}.$
(18)

We omit depth positions for simplicity. The first line of this equation (blue)
corresponds the AI group, while second and third lines (brown and red)
correspond to the CI group. Using the symmetricity of ${\bm{x}_{i}}$ and
${\bm{x}_{j}}$, the second summation part can be rewritten in a real form from
Eq. (18) as

$I(\bm{x^{\prime}})=\color[rgb]{0.2,0.2,0.6}\frac{w_{0}^{4}}{w^{4}}\sum_{i=1}^{N}B^{2}({\bm{x}_{i}})\exp{\left[{-\frac{4}{w^{2}}\left(\bm{x^{\prime}}-{\bm{x}_{i}}\right)^{2}}\right]}\\\
+\color[rgb]{0.43,0.21,0.1}\frac{2w_{0}^{4}}{w^{4}}\underset{(i>j)}{\sum_{i=1}^{N}\sum_{j=1}^{N}}B({\bm{x}_{i}})B({\bm{x}_{j}})\color[rgb]{0.8,0.0,0.0}\exp{\left[{-\frac{({\bm{x}_{i}}-{\bm{x}_{j}})^{2}}{w^{2}}}\right]}\exp{\left[{-\frac{4}{w^{2}}{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{i}}+{\bm{x}_{j}}}{2}}\right)}^{2}}\right]}\\\
\color[rgb]{0.8,0.0,0.0}\times\cos{\left[{(\phi_{i}-\phi_{j})+\frac{2nk_{0}}{R}{\left({{\bm{x}_{i}}-{\bm{x}_{j}}}\right)}\cdot{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{i}}+{\bm{x}_{j}}}{2}}\right)}}\right]}\color[rgb]{0,0,0}.\\\
$ (19)

By applying the same manipulations used to convert Eq. (6) to Eq. (10), which
we described in Section 3.3, we obtain the continuous form of the intensity-
OCT equation with the Gaussian PSF with defocus from Eq. (18) as

$I(\bm{x^{\prime}})=\color[rgb]{0.2,0.2,0.6}\frac{w_{0}^{4}}{w^{4}}\iint
D({\bm{x}_{p}})B^{2}({\bm{x}_{q}})\exp{\left[{-\frac{4}{w^{2}}\left(\bm{x^{\prime}}-{\bm{x}_{p}}\right)^{2}}\right]}d{\bm{x}_{p}}\\\
+\color[rgb]{0.43,0.21,0.1}\frac{w_{0}^{4}}{w^{4}}\iiiint
D^{\prime}({\bm{x}_{p}},{\bm{x}_{q}})B({\bm{x}_{p}})B({\bm{x}_{q}})\color[rgb]{0.8,0.0,0.0}\exp{\left[{-\frac{({\bm{x}_{p}}-{\bm{x}_{q}})^{2}}{w^{2}}}\right]}e^{i\phi_{pq}}\\\
\color[rgb]{0.8,0.0,0.0}\times\exp{\left[{-\frac{4}{w^{2}}{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{p}}+{\bm{x}_{q}}}{2}}\right)}^{2}}\right]}\exp
i{\left[{\frac{2nk_{0}}{R}{\left({{\bm{x}_{p}}-{\bm{x}_{q}}}\right)}\cdot{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{p}}+{\bm{x}_{q}}}{2}}\right)}}\right]}\color[rgb]{0.43,0.21,0.1}d{\bm{x}_{p}}d{\bm{x}_{q}}\\\
\color[rgb]{0.43,0.21,0.1}-\mathrm{the\,\,\,first\,\,\,line\,\,\,of\,\,\,R.H.S.},$
(20)

where R.H.S. means “the right hand side.” The last line, which is identical to
the AI group (blue), was subtracted from the second integral, so that the last
three lines of the equation (brown and red) represent the CI group. This
equation is the continuous-form imaging equation of an intensity-OCT image
with a complex Gaussian PSF with defocus. The first line and other lines
correspond to AI and CI, respectively.

Using the symmetricity of ${\bm{x}_{p}}$ and ${\bm{x}_{q}}$, and
$\phi_{pq}=-\phi_{qp}$, the second integral of Eq. (20) is rewritten in a real
form as

$I(\bm{x^{\prime}})=\color[rgb]{0.2,0.2,0.6}\frac{w_{0}^{4}}{w^{4}}\iint
D({\bm{x}_{p}})B^{2}({\bm{x}_{p}})\exp{\left[{-\frac{4}{w^{2}}\left(\bm{x^{\prime}}-{\bm{x}_{p}}\right)^{2}}\right]}d{\bm{x}_{p}}\\\
+\color[rgb]{0.43,0.21,0.1}\frac{w_{0}^{4}}{w^{4}}\iiiint
D^{\prime}({\bm{x}_{p}},{\bm{x}_{q}})B({\bm{x}_{p}})B({\bm{x}_{q}})\color[rgb]{0.8,0.0,0.0}\exp{\left[{-\frac{({\bm{x}_{p}}-{\bm{x}_{q}})^{2}}{w^{2}}}\right]}\exp{\left[{-\frac{4}{w^{2}}{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{p}}+{\bm{x}_{q}}}{2}}\right)}^{2}}\right]}\\\
\color[rgb]{0.8,0.0,0.0}\times\cos{\left[{\phi_{pq}+\frac{2nk_{0}}{R}{\left({{\bm{x}_{p}}-{\bm{x}_{q}}}\right)}\cdot{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{p}}+{\bm{x}_{q}}}{2}}\right)}}\right]}\color[rgb]{0.43,0.21,0.1}d{\bm{x}_{p}}d{\bm{x}_{q}}\\\
\color[rgb]{0.43,0.21,0.1}-\mathrm{the\,\,\,first\,\,\,line\,\,\,of\,\,\,R.H.S.}$
(21)

###  3.5 Interpretation of OCT intensity

Eqs. (6), (10), (19), and (21) provide an intuitive understanding of the
intensity-OCT image. All the equations suggest that the OCT intensity consists
of AI [i.e., the first summation or integral of each equation (blue) of Eqs.
(6), (10), (19), and (21)], and CI (i.e., the residual terms, brown and red).
AI forms a meaningful image, as we discuss in this paragraph, whereas CI
corresponds to speckle, which we discuss in the next paragraph. As can be seen
in the equations, for the intensity-OCT image formed by AI, the object to be
measured is $B^{2}({\bm{x}_{i}})$ for discrete forms and
$D(\bm{x})B^{2}(\bm{x})$ for continuous forms. Additionally, it is noteworthy
that this intensity-OCT imaging explicitly has a real- and positive-valued
PSF: ${P_{a}}^{2}$. We found that the AIs of the continuous-form equations
[Eqs. (10) and (21)] are the convolution of the real-positive PSF and the
object to be measured. And hence, the AI is equivalent to the incoherent image
of the object.

In contrast to AI, CI does not have an explicit PSF. To intuitively understand
the image property of CI, we split the product parts inside the summation or
integral into two parts: the object to be measured ($B^{2}({\bm{x}_{i}})$ or
$D(\bm{x})B^{2}(\bm{x})$) and the residual part (written in red). We refer to
the residual part as “pseudo-PSF” because it can be considered as the PSF for
each scatterer pair. The pseudo-PSF varies for each combination of scatterers;
hence, CI does not form a meaningful image. Namely, CI represents speckle.
Specifically for the intensity-OCT imaging under a Gaussian complex PSF [Eqs.
(19) and (21)], the pseudo-PSF is found to have the following characteristics.
First, the pseudo-PSF is the product of the Gaussian envelope (the second
exponential parts in CI) and a cosine wave. Second, the Gaussian envelope has
the same width as the explicit PSF of AI (the exponential part of AI). It is
consistent with a well-known fact; speckle size is a function of the OCT
resolution [5, 18]. Third, the center of the Gaussian is located at the center
of the two scatterers of each scatterer pair; that is, the center of the
Gaussian envelope walks around for each scatterer pair and is practically
unpredictable and random. Fourth, the frequency and phase offset, that is, the
shift of the cosine wave, depends on the positions of the scatterers of the
scatterer pair; hence, they are practically unpredictable and random.

These formulations indicate that the intensity-OCT image is a summation of a
speckle-free meaningful image, i.e., the incoherent image of the object, which
corresponds to AI, and a random speckle image corresponding to CI.

##  4 Volumetric differential contrast (VDC) imaging

###  4.1 Theory

In this section, we present a new OCT-based imaging method called volumetric
differential contrast (VDC) imaging, which we design using the DSM and a
similar theoretical framework presented in Section 3. This is an image
processing method that manipulates the en face OCT signal and can provide a
spatially differentiated image of a sample that highlights the sample
structure as contrast, similar to DIC. In contrast to standard DIC, VDC is
volumetric, that is, it provides a spatially differentiated image at each
depth of a sample.

An en face VDC image is obtained by multiplying an en face complex OCT by its
digitally shifted complex conjugate and then extracting the imaginary part. We
describe the formulation of VDC in this section.

In the first step of an en face VDC computation, we compute the product of an
en face complex OCT signal and the digitally laterally shifted complex
conjugate of this signal. Using the same theoretical flow that we used to
obtain Eq. (10), this product, which is an en face complex signal, is written
as

$s{\left({\bm{x^{\prime}}+\frac{\Delta\bm{x}}{2}}\right)}s^{*}{\left({\bm{x^{\prime}}-\frac{\Delta\bm{x}}{2}}\right)}=\color[rgb]{0.2,0.2,0.6}\iint
D({\bm{x}_{p}})B^{2}({\bm{x}_{p}}){{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{p}}+\frac{\Delta\bm{x}}{2}\right)}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{p}}-\frac{\Delta\bm{x}}{2}\right)}\\\
\color[rgb]{0.2,0.2,0.6}\times\exp{i\left[\varphi\left(\bm{x^{\prime}}-{\bm{x}_{p}}+\frac{\Delta\bm{x}}{2}\right)-\varphi\left(\bm{x^{\prime}}-{\bm{x}_{p}}-\frac{\Delta\bm{x}}{2}\right)\right]}d{\bm{x}_{p}}\\\
+\color[rgb]{0.43,0.21,0.1}\iiiint
D^{\prime}({\bm{x}_{p}},{\bm{x}_{q}})B({\bm{x}_{p}})B({\bm{x}_{q}})\color[rgb]{0.8,0.0,0.0}e^{i(\phi_{pq})}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{p}}+\frac{\Delta\bm{x}}{2}\right)}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{q}}-\frac{\Delta\bm{x}}{2}\right)}\\\
\qquad\color[rgb]{0.8,0.0,0.0}\times\exp{i\left[\varphi\left(\bm{x^{\prime}}-{\bm{x}_{p}}+\frac{\Delta\bm{x}}{2}\right)-\varphi\left(\bm{x^{\prime}}-{\bm{x}_{q}}-\frac{\Delta\bm{x}}{2}\right)\right]}\color[rgb]{0.43,0.21,0.1}d{\bm{x}_{p}}d{\bm{x}_{q}}\\\
\color[rgb]{0.43,0.21,0.1}-\mathrm{the\,\,\,first\,\,\,integral\,\,\,of\,\,\,R.H.S.}\color[rgb]{0,0,0},$
(22)

where $\Delta\bm{x}$ is the digital shift amount between the two signals,
typically one or two lateral pixels, and should be smaller than the spot size
of the probe beam. The letter colors corresponds to those in the previous
equations. The detailed derivation of this equation is in the Appendix.

By assuming the Gaussian complex PSF described in Section 3.4, that is, by
substituting Eqs. (16) and (17) into Eq. (22), we obtain the en face complex
signal:

$s{\left({\bm{x^{\prime}}+\frac{\Delta\bm{x}}{2}}\right)}s^{*}{\left({\bm{x^{\prime}}-\frac{\Delta\bm{x}}{2}}\right)}=\color[rgb]{0.2,0.2,0.6}\frac{w_{0}^{4}}{w^{4}}\exp{\left({-\frac{\Delta\bm{x}^{2}}{w^{2}}}\right)}\iint
D({\bm{x}_{p}})B^{2}({\bm{x}_{p}})\\\
\color[rgb]{0.2,0.2,0.6}\times\exp{\left[{-\frac{4}{w^{2}}(\bm{x^{\prime}}-{\bm{x}_{p}})^{2}}\right]}\exp
i{\left[{-\frac{2nk_{0}}{R}\Delta\bm{x}\cdot{\left({\bm{x^{\prime}}-{\bm{x}_{p}}}\right)}}\right]}d{\bm{x}_{p}}\\\
\quad+\color[rgb]{0.43,0.21,0.1}\frac{w_{0}^{4}}{w^{4}}\exp{\left({-\frac{\Delta\bm{x}^{2}}{w^{2}}}\right)}\iiiint
D^{\prime}({\bm{x}_{p}},{\bm{x}_{q}})B({\bm{x}_{p}})B({\bm{x}_{q}})\\\
\color[rgb]{0.8,0.0,0.0}\times\exp{\left[{-\frac{2}{w^{2}}({\bm{x}_{p}}-{\bm{x}_{q}})\cdot{\left({\frac{{\bm{x}_{p}}+{\bm{x}_{q}}}{2}-\Delta\bm{x}}\right)}}\right]}e^{i\phi_{pq}}\\\
\color[rgb]{0.8,0.0,0.0}\times\exp{\left[{-\frac{4}{w^{2}}{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{p}}+{\bm{x}_{q}}}{2}}\right)}^{2}}\right]}\exp
i{\left[{\frac{2nk_{0}}{R}{\left({{\bm{x}_{p}}-{\bm{x}_{q}}-\Delta\bm{x}}\right)}\cdot{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{p}}+{\bm{x}_{q}}}{2}}\right)}}\right]}\color[rgb]{0.43,0.21,0.1}d{\bm{x}_{p}}d{\bm{x}_{q}}\\\
\color[rgb]{0.43,0.21,0.1}-\mathrm{the\,\,\,first\,\,\,integral\,\,\,of\,\,\,R.H.S.}$
(23)

Similar to the conventional intensity-OCT image [Eq. (10)], this en face
signal consists of the AI part [the first two lines of Eq. (23), blue] and the
CI part (the next four lines, brown and red). The AI part forms an imaging
formula in which the object to be measured is the same as the conventional
intensity-OCT, that is, $D(\bm{x})B^{2}(\bm{x})$, the product of the scatterer
density and the squared scattering potential. The PSF of this imaging formula
is a complex PSF (the second line of the equation), which is the product of a
Gaussian and a lateral linear phase slope. By contrast, the CI part is
practically random, which is, again, similar to the conventional intensity-
OCT.

Now we define the VDC image as the imaginary part of Eq. (23) as

$\text{Im}\left[s{\left({\bm{x^{\prime}}+\frac{\Delta\bm{x}}{2}}\right)}s^{*}{\left({\bm{x^{\prime}}-\frac{\Delta\bm{x}}{2}}\right)}\right]=\color[rgb]{0.2,0.2,0.6}\frac{w_{0}^{4}}{w^{4}}\exp{\left({-\frac{\Delta\bm{x}^{2}}{w^{2}}}\right)}\iint
D({\bm{x}_{p}})B^{2}({\bm{x}_{p}})\\\
\color[rgb]{0.2,0.2,0.6}\times\exp{\left[{-\frac{4}{w^{2}}(\bm{x^{\prime}}-{\bm{x}_{p}})^{2}}\right]}\sin{\left[{-\frac{2nk_{0}}{R}\Delta\bm{x}\cdot{\left({\bm{x^{\prime}}-{\bm{x}_{p}}}\right)}}\right]}d{\bm{x}_{p}}\\\
\quad+\color[rgb]{0.43,0.21,0.1}\frac{w_{0}^{4}}{w^{4}}\exp{\left({-\frac{\Delta\bm{x}^{2}}{w^{2}}}\right)}\iiiint
D^{\prime}({\bm{x}_{p}},{\bm{x}_{q}})B({\bm{x}_{p}})B({\bm{x}_{q}})\\\
\color[rgb]{0.8,0.0,0.0}\times\exp{\left[{-\frac{2}{w^{2}}({\bm{x}_{p}}-{\bm{x}_{q}})\cdot{\left({\frac{{\bm{x}_{p}}+{\bm{x}_{q}}}{2}-\Delta\bm{x}}\right)}}\right]}\exp{\left[{-\frac{4}{w^{2}}{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{p}}+{\bm{x}_{q}}}{2}}\right)}^{2}}\right]}\\\
\color[rgb]{0.8,0.0,0.0}\times\sin{\left[{\phi_{pq}+\frac{2nk_{0}}{R}{\left({{\bm{x}_{p}}-{\bm{x}_{q}}-\Delta\bm{x}}\right)}\cdot{\left({\bm{x^{\prime}}-\frac{{\bm{x}_{p}}+{\bm{x}_{q}}}{2}}\right)}}\right]}\color[rgb]{0.43,0.21,0.1}d{\bm{x}_{p}}d{\bm{x}_{q}}\\\
\color[rgb]{0.43,0.21,0.1}-\mathrm{the\,\,\,first\,\,\,integral\,\,\,of\,\,\,R.H.S.}$
(24)

In this equation, the phase parts of Eq. (23), i.e., the last exponential part
of each integral, became sinusoidal function. The last four lines of Eq. (24)
(brown and red) correspond to the CI part. Similar to the case of the
intensity-OCT formula, $\phi_{pq}$ and pseudo-PSF (red) are practically
random; hence, CI becomes a random speckle pattern. By omitting this random CI
part, we can define the VDC image
[$\mathrm{VDC}{\left({\bm{x^{\prime}}}\right)}$] and write it as the first
integral of Eq. (24):

$\mathrm{VDC}{\left({\bm{x^{\prime}}}\right)}\equiv\color[rgb]{0.2,0.2,0.6}\frac{w_{0}^{4}}{w^{4}}\exp{\left({-\frac{\Delta\bm{x}^{2}}{w^{2}}}\right)}\\\
\color[rgb]{0.2,0.2,0.6}\times\iint
D(\bm{x})B^{2}(\bm{x})\exp{\left[{-\frac{4}{w^{2}}{\left({\bm{x^{\prime}}-\bm{x}}\right)}^{2}}\right]}\sin\left[-\frac{2nk_{0}}{R}\Delta\bm{x}\cdot{\left({\bm{x^{\prime}}-\bm{x}}\right)}\right]d\bm{x}\color[rgb]{0,0,0}.$
(25)

We interpret this imaging formula as the convolution of the object
$D(\bm{x})B^{2}(\bm{x})$ with an imaging PSF that is a Gaussian-windowed sine
function. This imaging PSF function is similar to the response function (i.e.,
the imaging PSF) of DIC [19] and also a differential operator used in digital
image processing [20]. Hence, it works as a differential operator in VDC
imaging and provides the differential contrast of the sample structure.

Figure 3:  Examples of the imaging PSFs of VDC imaging with different shift
amounts (a) and a different defocus amount (b). The one-dimensional plots are
the profile of the two-dimensional en face PSF at the center of the PSF in the
shift direction. (a) The larger the shift amount, the higher the frequency of
the sine part of the PSF and vice versa. When the shift is zero, the sine part
becomes zero; hence, the image is not formed. The shift amount does not affect
the Gaussian part. (b) The larger the defocus amount, the wider the Gaussian
width and the higher the frequency of the sine part, and vice versa, provided
that the defocus amount is larger than the Rayleigh length. When the defocus
amount is zero, the sine part becomes zero again; hence, the image is not
formed. The appropriate selection of the shifting and defocus amounts leads to
the PSF becoming a differential operator. The direction of differentiation can
be flipped by flipping the shift and defocus directions. Note that the imaging
PSF of VDC is analogous to the response function of DIC, and the shear amount
and direction of DIC are analogous to the size and direction of the imaging
PSF of VDC.

The shape of the imaging PSF of VDC (i.e., the differential operator) is
altered by the defocus amount $z-z_{0}$ and shift $\Delta\bm{x}$. Figure 3
exemplifies several PSFs with different shifts (a) and defocus amounts (b). We
computed PSFs using OCT system parameters identical to those in the
experimental setup described in Section 4.2: $w_{0}=2.41$ $\muup\text{m}$,
$n=1.38$, and $k_{0}=7.48\times 10^{6}$ rad/m. The shifting amounts were the
integer multiples of the lateral pixel size (1.95 $\muup\text{m}$) and the
defocus amounts were expressed as the distance from the focal depth, that is,
$z-z_{0}$ in Eqs. (13) and (14). For Fig. 3(a), we set the defocus amount to
200 $\muup\text{m}$, whereas for Fig. 3(b), we set the shifting amount to 1.95
$\muup\text{m}$. The one-dimensional plots are the profiles at the center of
the PSF, i.e., $\bm{x^{\prime}}=\bm{x}$, and are in the digital shifting
direction. The columns “Gaussian part” and “sine part” in the figure represent
the Gaussian and sine parts within the integral of Eq. (25), respectively.

As shown in Fig. 3(a), the shift affected the frequency of the sine function
in the imaging PSF and did not affect the Gaussian function. By contrast, the
defocus amount affected both the Gaussian and sine functions of the imaging
PSF [Fig. 3(b)]. Please note that the absolute defocus amount was larger than
the Rayleigh length (21.7 $\muup\text{m}$) in the examples in the figure.

By properly selecting the shift and defocus amounts, we can obtain a good
differential operator with an arbitrary differential distance, which is
defined as the distance between two main positive and negative peaks. The
differential distances are shown at the top right of the 2D-PSF images. (By
contrast, an example of the “inappropriate” selection of these amounts is a
shift that is too large, which fluctuates the operator within the Gaussian and
an appropriate differential contrast cannot be obtained.) When the defocus
direction is reversed, that is, the plus/minus sign of the defocus is changed,
the sine function and differential direction are flipped, as shown in Fig.
3(b). When the shift is zero, the sine function becomes zero; hence,
$\mathrm{VDC}(\bm{x^{\prime}})=0$; that is, differential contrast is not
obtained. Similarly, when the defocus amount is zero, that is, $z=z_{0}$, $R$
becomes infinity and the sine part of the operator becomes $\sin 0=0$; hence,
contrast is not obtained again. Thus, not only the shift but also the defocus
is mandatory for VDC.

In our particular implementation, we computationally apply defocusing as
detailed in Section 4.3.3; hence, both the shift and defocus amount can be
selected after image acquisition. It is also noteworthy that the size and
direction of the differential operator are analogous to the shear amount and
the shear direction of DIC imaging. Because the shift and defocus amount can
be set after image acquisition, VDC can generate multiple DIC-like images with
arbitrary shear amounts and directions using post-acquisition signal
processing.

###  4.2 Experimental validation

####  4.2.1 System, samples, and image acquisition protocol

For the experimental validation of VDC imaging, we used a fiber-based
spectral-domain OCT system[14]. The light source was a superluminescent diode
(SLD, Superlum, Ireland) with a center wavelength and bandwidth of 840 nm and
100 nm, respectively. In the sample arm, the probe beam with a diameter of 4.0
mm passed through an objective (LSM02BB, Thorlabs, NJ) with a focal length of
18 mm and incident on the sample. The axial resolution was 3.8 $\muup\text{m}$
in tissue and the in-focus lateral resolution was 4.9 $\muup\text{m}$. The
spectral interference signal was acquired using a spectrometer (Cobra S800,
Wasatch Photonics, NC). The scanning speed was 50,000 lines/s with 2,048
spectral pixels.

To reconstruct the OCT volume, spectral rescaling, numerical dispersion
compensation, and fixed-pattern noise removal were applied. The depth-
independent bulk phase error caused by system imperfection and air turbulence
was numerically corrected after the reconstruction of the OCT volume using the
smart-integration-path method [21]. Note that, since our image processing
exploits complex signals, high phase stability is required

Figure 4:  Cross-sectional (a, c) and en face (b, d) intensity-OCT images of a
zebrafish (a, b) and tumor spheroid (c, d). The yellow lines in the cross-
sectional images indicate the depth positions of the en face images.
Similarly, the yellow lines in the en face images indicate the positions of
the cross-sectional images.

An 8-day-old zebrafish larva and a human breast cancer (MCF7) spheroid were
measured. The zebrafish was anesthetized and embedded in agarose gel. The
lateral scanning areas were 1 mm $\times$ 1 mm for the zebrafish and 0.8 mm
$\times$ 0.8 mm for the spheroid, with 512 $\times$ 512 A-lines. Figure 4
shows examples of the cross-sectional and en face OCT images of the (a)
zebrafish and (b) spheroid. The yellow dotted lines in each image indicate the
depths of the en face images and position of the cross-sectional images. These
datasets were used for the validations in Section 4.3.

####  4.2.2 Processing of VDC images

Figure 5:  Processing flow of VDC imaging. (1) A complex en face slice was
extracted from a complex OCT volume. (2) Computational refocusing and
defocusing were applied to obtain an arbitrarily defocused signal ($s_{1}$).
(3) A digitally shifted copy of the en face signal was generated ($s_{2}$) and
(4) complex conjugated. (5) The product of $s_{1}$ and the complex conjugate
of $s_{2}$ was computed. (6) VDC was obtained as the imaginary part of the
product ($s_{1}s_{2}^{*}$).

The VDC processing flow is summarized in Fig. 5. (1) An en face plane was
extracted from a complex OCT volume. (2) The en face signal was first
refocused using computational refocusing based on a forward-propagation model
[22, 21] and then arbitrary numerical defocusing was applied (similar to that
in Section 5.1.1 of Ref. [23]). (3) This complex en face OCT signal with
arbitrary defocusing, $s_{1}(\bm{x^{\prime}})$, was copied and digitally
shifted in the lateral direction to obtain $s_{2}(\bm{x^{\prime}})$. (4, 5)
Then the product of $s_{1}(\bm{x^{\prime}})$ and the complex conjugate of
$s_{2}(\bm{x^{\prime}})$ was computed. (6) Finally, the en face VDC image was
obtained as the imaginary part of the product.

###  4.3 Results

####  4.3.1 Shifting direction dependency of VDC imaging

Figure 6:  En face VDC images of the (a) zebrafish and (b) spheroid with
multiple shift directions. The edge structures normal to the differential
direction (i.e., shift direction) are emphasized. The insets show
corresponding imaging PSFs. The plots on the right show some representative
line profiles of the VDC images. The positions of the line profiles are
indicated by long yellow arrows on the VDC images.

The En face VDC images of the (a) zebrafish and (b) spheroid are summarized in
Fig. 6, which demonstrates the multi-directional differential imaging
capability of VDC. The plots on the right show some representative line
profiles of the VDC images. The positions of the line profiles are indicated
by long yellow arrows on the VDC images. The shifting amounts
$\Delta\bm{x}=(\Delta x,\Delta y)$ were the combination of -1, 0, and +1
pixels (i.e., -1.95, 0.0, and +1.95 $\muup\text{m}$ for the zebrafish and
-1.56, 0.0, and +1.56 $\muup\text{m}$ for the spheroid) for the $x$ and $y$
directions. These shift amounts were selected to be the integer multiples of
the pixel size to avoid sub-pixel interpolation of the complex OCT image. The
defocus amount was +200 $\muup\text{m}$ for all images. Note that this defocus
was not a physical defocus but a numerically induced defocus by using complex
signal processing as mentioned in Section 4.2.2. The combination of the shift
and defocus amounts were selected to make the imaging PSF as a good
differential operator. The small insets represent the imaging PSF.

As predicted by the imaging PSF shape of VDC, these images looked like the DIC
image. In the VDC images, the edge structures normal to the differential
direction (shift direction) were particularly emphasized. As expected, the VDC
image was not obtained when there was no shift. Note that, because DIC cannot
be performed on a thick sample, such as the zebrafish and spheroid, a direct
comparison of VIC and DIC cannot be performed.

####  4.3.2 Shifting amount dependency of VDC imaging

Figure 7:  En face VDC images of the (a)–(e) zebrafish and (f)–(j) spheroid
with shifting amounts of -2, -1, 0, +1, and +2 pixels. The small insets show
the corresponding imaging PSFs and the values above them show the differential
distances. $\pm$1-pix shift yielded a larger differential distance than
$\pm$2-pix shift, and hence emphasized larger structures than $\pm$2-pixel-
shift images.

By altering the shifting amount, VDC images with multiple differential
distances are obtained as shown in Fig. 7. The shifting amounts were -2, -1,
0, +1, and +2 pixels (i.e., -3.9, -1.95, 0.0, +1.95, and +3.9 $\muup\text{m}$
for the zebrafish and -3.12, -1.56, 0.0, +1.56, +3.12 $\muup\text{m}$ for the
spheroid) along the $x$-axis. The defocus amount was +200 $\muup\text{m}$ for
all images. The insets show the corresponding imaging PSFs and differential
distances are shown above.

$\pm$1-pixel shifts [Fig. 7(b), (d), (g), and (i)] yielded larger differential
distances than $\pm$2-pixel shifts [Fig. 7(a), (e), (f), and (j)]. Hence,
$\pm$1-pixel-shift images emphasized larger structures than $\pm$2-pixel-shift
images. It is also noteworthy that the imaging PSFs of $\pm$2-pixel shifts
exhibited small side peaks (ripples). This may have deteriorated the VDC
images. When the shifting amount was zero, contrast was not obtained as
expected [Fig. 7(c) and (h)].

####  4.3.3 Defocus amount dependency of VDC imaging

Figure 8:  En face VDC images of the (a)–(e) zebrafish and (f)–(j) spheroid
with defocus amounts of -200, -100, 0, +100, and +200 $\muup\text{m}$. The
small insets show the corresponding imaging PSFs and the values above them
show the differential distances. The differential distance increased as the
amount of defocus increased. The differential direction reversed when the
defocus direction reversed.

The effect of the defocus amount is demonstrated in Fig. 8, where the defocus
amounts were -200, -100, 0, +100, and +200 $\muup\text{m}$. The shifting
amount was (+1, 0) pixel for the $(x,y)$-direction for all images. The first
and second rows show the VDC images of the zebrafish and spheroid,
respectively. The insets show the corresponding imaging PSFs and differential
distances are shown above.

Increasing the absolute defocus amount increased the differential distance.
When the defocusing direction was reversed, the differential direction was
also reversed. The VDC contrast became very low for the in-focus
(0-$\muup\text{m}$ defocus) and [Fig. 8(c) and (e)]. It should be noted that,
in theory, the VDC is expected to vanish at in-focus, but low contrast was
visible in the images. This residual contrast might be caused by the
imperfection of numerical refocusing.

####  4.3.4 Volumetric differential imaging

Figure 9:  Volumetrically processed VDC images of the zebrafish and tumor
spheroid. (a1–f1) En face VDC images at each depth of the zebrafish. The depth
positions are indicated in the cross-sectional OCT image (g1). (a2–g2) Same
images for the tumor spheroid. The fly-through en face images are available as
Visualization-1 (zebrafish) and Visualization-2 (tumor spheroid).

As shown above, VDC yielded DIC-like differential images. However, unlike DIC,
VDC achieved volumetric differential imaging by processing each en face plane
in an OCT volume. Figure 9 shows examples of the zebrafish (a1–g1) and tumor
spheroid (a2–g2). (a1)–(f1) and (a2)–(f2) are the en face VDC images extracted
from the VDC volumes (from the upper side to the bottom side). The depth
locations of the en face images are indicated in the OCT cross-sections (g1)
and (g2). The defocus and shifting amounts were +200 $\muup\text{m}$ and +1
pix along the $x$-axis for all images. The fly-through en face images for all
depths are available as Visualization-1 (zebrafish) and Visualization-2 (tumor
spheroid).

##  5 Discussion

###  5.1 Advantage of VDC imaging

As shown in Section 4, VDC imaging provides images similar to DIC, but VDC
imaging has two advantages over DIC. First, VDC is volumetric and can image
thick tissue. Second, VDC can generate images with several differential
conditions with a single acquisition; that is, the differential direction and
differential distance, which correspond to the shear direction and amount of
DIC, respectively, can be controlled by the parameters of post-acquisition
signal processing. This means that VDC can be regarded as a volumetric DIC
with post-hoc control of the shear direction and amount.

###  5.2 Difference between VDC and DIC

Despite their image similarity, the VDC and DIC have some fundamental
differences in their contrast mechanisms. DIC’s contrast is mainly from the
optical path length difference between two adjacent lateral points. In
addition, DIC is a transillumination microscope, and the sample is thin, low
scattering, and typically, transparent. And hence, the DIC’s contrast is
mainly based on the refractive index distribution of the sample.

On the other hand, VDC is obtained by OCT, which is a reflection (back-
scattering) measurement modality, and the sample is assumed to exhibit high
scattering. And hence, VDC’s contrast might be mainly based on the back-
scattering.

One of the not-well-described effects in the current VDC theory is the effect
of superior tissues. VDC may be sensitive to the refractive index of the
superior tissues, because it alters the phase of the complex en face OCT
signal at the depth of interest. However, the present theory only considers an
image at a particular depth, and does not formulate the cumulative effect of
the superior tissues. Extension of VDC theory to account for the cumulative
effect is one of the future work.

In addition, three-dimensional imaging PSF of VDC has not yet been
investigated. As demonstrated by Villiger et al.[24] and Ralston et al.[25],
the three-dimensional imaging property of OCT can be formulated. The extension
of VDC imaging theory into three-dimensional space is also one of the future
works.

Well-characterized-phantom based experimental validation might be helpful to
further understand the imaging property of VDC. However, for this validation,
an inhomogeneous scattering phantom that has a slowly varying refractive index
distribution and multiple domains with different scatterer densities are
necessary. The establishment of such characterized-phantoms and further
experimental validation of VDC might be one of the future studies.

###  5.3 Intensity based differential image generation

Figure 10: Differential images generated from intensity based signal
processing. The images were generated by convolving the linear-scale OCT
intensity image with the imaging PSF of Fig. 6 which is shown also in the
insets of this figure. The images corresponds to those of Fig. 6. Since this
operation differentiate the meaningful OCT image and speckle by the same
differentiation operator, the signals originated from the meaningful OCT image
and the speckle are hardly distinguished.

It is noteworthy that, as far as we consider only the AI part, the
differential contrast image can be obtained by one of the standard intensity-
based image processing methods. Figure 10 shows the differential contrast
images obtained by the intensity-based image processing corresponding to Fig.
6. Here the images were obtained by convolving the linear-scale OCT intensity
image with the imaging PSF of Fig. 6, and we denote these images as “intensity
differential images.”

Two differences are noted between VDC (Fig. 6) and the intensity differential
images (Fig. 10). First, the resolution of this intensity differential image
might be slightly worse than that of VDC. It is because the original
intensity-OCT image used for the intensity differential image has already been
affected by the intensity PSF [$P_{a}^{2}(\bm{x})$ of Eq. (10)]. However, the
imaging PSF of VDC is much larger than the intensity PSF due to the numerical
defocusing. The defocus amount was 200 $\muup\text{m}$ in this particular
case, while the Rayleigh length is 21.7 $\muup\text{m}$. And hence, this
reduction of the resolution might be negligible.

The second and more significant difference is that the intensity differential
images exhibit more granular appearance than the VDC images. This difference
can be explained by the effect of the CI part. In the intensity differential
images, CI part, i.e., speckle is also differentiated by the same
differentiation kernel with the meaningful AI part. And hence, “differential
speckle image” is overlaid on the meaningful differential images. Since both
of them have the same resolution, it is hard to distinguish the differential
speckle from the meaningful differential image.

On the other hand, the CI part of the VDC [i.e., the last two lines of Eq.
(24)] does not have an evident PSF. It has only a pseudo-PSF [see the
paragraph just after Eq. (24)], and the width and the magnitude of pseudo-PSF
alters for all combinations of the contributing scatterers. And hence, the
speckle-originated signal (CI part) appears quite differently from the
meaningful differential image (AI part) in VDC. Quantitative comparison of
visibility of specific image features in specific applications between VDC
intensity-based image-processing methods might be one of the future works in
the transnational research.

###  5.4 Artifacts in Doppler OCT

Doppler OCT enables flow imaging in tissue[26]. Among several variations of
Doppler OCT, conventional versions compute the phase difference between
adjacent A-lines in a B-scan to obtain a Doppler OCT image [27, 28]. This
version of Doppler OCT is known to exhibit a phase noise pattern caused by the
structure of the sample, that is, the structural phase artifact. This
structural phase artifact can be mathematically described by extending the
formulation of VDC as follows:

This version of Doppler OCT can be written by modifying Eq. (23) to
$s{\left({\bm{x^{\prime}}+\Delta\bm{x}/2}\right)}$
$s^{*}{\left({\bm{x^{\prime}}-\Delta\bm{x}/2}\right)}$ $e^{i\varphi_{D}}$,
where the newly added phase $\varphi_{D}$ is the Doppler phase shift and
$\Delta\bm{x}$ is the separation of the adjacent A-lines. The AI part of this
modified equation becomes

$\mathrm{AI}=\frac{w_{0}^{4}}{w^{4}}\exp{\left({-\frac{\Delta\bm{x}^{2}}{w^{2}}}\right)}\exp{\left({i\varphi_{D}}\right)}\iint
D({\bm{x}_{p}})B^{2}({\bm{x}_{p}})\\\
\times\exp{\left[{-\frac{4}{w^{2}}(\bm{x^{\prime}}-{\bm{x}_{p}})^{2}}\right]}\exp
i{\left[{-\frac{2nk_{0}}{R}\Delta\bm{x}\cdot{\left({\bm{x^{\prime}}-{\bm{x}_{p}}}\right)}}\right]}d{\bm{x}_{p}}$
(26)

and the Doppler OCT signal is defined as the phase of this equation.

This formulation raises two points about the structural phase artifact. First,
if the static structures of the adjacent A-lines are identical, i.e.,
$s(\bm{x^{\prime}}+\Delta\bm{x}/2)=s(\bm{x^{\prime}}-\Delta\bm{x}/2)$, the
Doppler OCT signal becomes the Doppler phase shift as

$\mathrm{angle}\left[s{\left({\bm{x^{\prime}}+\frac{\Delta\bm{x}}{2}}\right)}s^{*}{\left({\bm{x^{\prime}}-\frac{\Delta\bm{x}}{2}}\right)}e^{\varphi_{D}}\right]=\varphi_{D}.$
(27)

Namely, no structural phase artifact appears if there is no structure.

Second, if there is no defocus, i.e., $R\rightarrow\infty$ in Eq. (26), the
phase of Eq. (26) (the Doppler OCT signal) becomes the Doppler phase;
otherwise, the Doppler OCT signal is not identical to the Doppler phase. It is
noteworthy that, for the second point, we did not assume
$s(\bm{x^{\prime}}+\Delta\bm{x}/2)=s(\bm{x^{\prime}}-\Delta\bm{x}/2)$, i.e.,
the structural phase artifact does not appear, even if there is a structure as
far as there is no defocus. To summarize, the structural phase artifact
appears only when the structure and defocus coexist.

Another structure-related issue is known for Doppler OCT; that is, “structural
decorrelation,” which means that the measurement accuracy of the Doppler phase
degrades as the distance between the adjacent A-lines ($\Delta\bm{x}$)
increases. This effect can be related to the first exponential part of Eq.
(26); that is, the AI signal decreases as $\Delta\bm{x}$ increases and leads
to a low signal-to-noise ratio, and hence, low phase measurement accuracy.

###  5.5 Defocusing effect of computational multi-directional (CMD)-OCT

Oida et al.demonstrated point-scanning computational multi-directional
(CMD)-OCT for qualitatively visualizing the microscopic phase structure of the
sample [14]. This method exploits the phase angle of
$s{\left({\bm{x^{\prime}}+\Delta\bm{x}/2}\right)}s^{*}{\left({\bm{x^{\prime}}-\Delta\bm{x}/2}\right)}$,
and hence is similar to our VDC method.

Although it was empirically known that the contrast of CMD-OCT becomes very
low if the image is in focus [15], its theoretical basis was not known. This
low contrast has now been explained by Eq. (23) of the VDC formulation; that
is, if there is no defocus (i.e., $R\rightarrow\infty$), the phase of the AI
part (meaningful part) of
$s{\left({\bm{x^{\prime}}+\Delta\bm{x}/2}\right)}s^{*}{\left({\bm{x^{\prime}}-\Delta\bm{x}/2}\right)}$
becomes zero. This fact explains the contrast vanishing of CMD-OCT at in-
focus.

###  5.6 The relation of our and previous OCT imaging theories

Zhou et al.has formulated the OCT signal as the collection of contributions
from discrete scatterers (Section 7 of [12]), and gave a similar
representation with our discrete equation of intensity-OCT [Eq. (6)]. Namely,
they represented the OCT image and speckle as two independent additive
components.

Our work added two points to their work. At first, we derived a continuous
form of the intensity-OCT imaging equation [Eq. (10)]. This equation clarified
that the meaningful OCT image is the convolution of the object
[$D(\bm{x})B^{2}(\bm{x})$] and a PSF of OCT intensity image
$P^{2}_{a}(\bm{x})$. Second, our work accounts for the phase of complex PSF.
It has enabled the design of VDC (Section 4), interpretations of artifacts in
Doppler OCT (Section 5.4), and understand the imaging property of
computational CMD-OCT (Section 5.5).

Our OCT formulation does not represent the three-dimensional imaging property,
similar to that discussed in Section 5.2 for VDC. It is because we have
assumed that the lateral and axial imaging properties are independent of each
other. Namely, low numerical aperture (NA) was tacitly assumed. Villiger et
al.[24] and Ralston et al.[25] have formulated the three-dimensional imaging
property of OCT, in which the lateral and axial imaging properties interact
with each other. The extension of our OCT formulation to account for the axial
imaging property and its interaction with the lateral imaging property might
be a future work.

Another point to be discussed is the applicability of our formulation strategy
to full-field (FF-) OCT. There are several works of formulation of FF-OCT
imaging. For example, Tricoli et al.formulated imaging property of partially
coherent FF-OCT [29], and Barolle et al.gave a theoretical model which
recapitulates the manifestation of aberrations in FF-OCT [30]. Our formulation
is potentially adoptable to FF-OCT. However, after Eq. (15), it is assumed
that the illumination and collection pupil functions are identical. And this
assumption is valid for point-scanning OCT but not true for FF-OCT with
spatially incoherent light sources. We can make our theory applicable to FF-
OCT by reformulating it from Eq. (10).

##  6 Conclusion

We presented a new mathematical model called DSM to represent a sample to be
measured by OCT. Using this model, we reformulated lateral OCT imaging and
also designed a new imaging method called VDC imaging.

Using the OCT reformulation, the meaningful signal and speckle were
successfully formulated as two different mathematical entities. We also found
that the meaningful signal and speckle originated from the auto- and cross-
numerical interactions of the scatters’ signal contribution, respectively.

VDC imaging is a combination of physical OCT imaging and subsequent complex
signal processing, and we derived an imaging formula that expresses the
physical and numerical process simultaneously. The imaging PSF that appeared
in this formula suggested that VDC imaging yields an image similar to that
yielded by DIC. We validated the VDC method by measuring biological samples.

Finally, the formulation strategy has the newly explained mathematical
mechanism of the well-known but previously not well-explained artifacts of
Doppler OCT, and the imaging property of CMD-OCT.

The presented formulation strategy may be applicable to other OCT-based
imaging methods in the future. Although the present formulation focuses on the
lateral imaging property, a similar approach based on DSM can be used to
extend this formulation for axial and three-dimensional imaging properties.

## Funding

Core Research for Evolutional Science and Technology (JPMJCR2105); Japan
Society for the Promotion of Science (21H01836, 22K04962, 22F22355, 22H05556,
21H02997); Japan Science and Technology Agency (JPMJMI18G8); Austrian Science
Fund (Schrödinger grand J4460).

## Disclosures

Tomita, Makita, Morishita, Abd El-Sadek, Mukherjee, Lichtenegger, Yasuno:
Nikon (F), Sky Technology (F), Yokogawa Electric Corp. (F), Kao Corp. (F),
Topcon (F). Fukutake: Nikon (E). Tamaoki, Bian, Kobayashi: None. Mori,
Matsusaka: None.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Appendix

In this appendix, we derive Eq. (22) in Section 4. The derivation is based on
a logical flow similar to that used to obtain Eq. (10) from Eq. (5). First, we
calculate the product of an en face complex OCT signal and its digitally
laterally shifted complex conjugate as

$\displaystyle
s{\left({\bm{x^{\prime}}+\frac{\Delta\bm{x}}{2}}\right)}s^{*}{\left({\bm{x^{\prime}}-\frac{\Delta\bm{x}}{2}}\right)}$
$\displaystyle\qquad=\left[\sum_{i=1}^{N}B({\bm{x}_{i}})e^{i\phi_{i}}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{i}}+\frac{\Delta\bm{x}}{2}\right)}\,e^{i\varphi{\left({\bm{x^{\prime}}-{\bm{x}_{i}}+\frac{\Delta\bm{x}}{2}}\right)}}\right]$
$\displaystyle\qquad\quad\times\left[\sum_{j=1}^{N}B({\bm{x}_{j}})e^{-i\phi_{j}}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{j}}-\frac{\Delta\bm{x}}{2}\right)}\,e^{-i\varphi{\left({\bm{x^{\prime}}-{\bm{x}_{j}}-\frac{\Delta\bm{x}}{2}}\right)}}\right]$
$\displaystyle\qquad=\color[rgb]{0.2,0.2,0.6}\sum_{i=1}^{N}B^{2}({\bm{x}_{i}}){{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{i}}+\frac{\Delta\bm{x}}{2}\right)}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{j}}-\frac{\Delta\bm{x}}{2}\right)}e^{i\left[\varphi{\left({\bm{x^{\prime}}-{\bm{x}_{i}}+\frac{\Delta\bm{x}}{2}}\right)}-\varphi{\left({\bm{x^{\prime}}-{\bm{x}_{i}}-\frac{\Delta\bm{x}}{2}}\right)}\right]}$
$\displaystyle\qquad\quad+\color[rgb]{0.43,0.21,0.1}\underset{(i\neq
j)}{\sum_{i=1}^{N}\sum_{j=1}^{N}}B({\bm{x}_{i}})B({\bm{x}_{j}})e^{i(\phi_{i}-\phi_{j})}$
$\displaystyle\qquad\qquad\color[rgb]{0.43,0.21,0.1}\times\color[rgb]{0.8,0.0,0.0}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{i}}+\frac{\Delta\bm{x}}{2}\right)}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{j}}-\frac{\Delta\bm{x}}{2}\right)}e^{i\left[\varphi{\left({\bm{x^{\prime}}-{\bm{x}_{i}}+\frac{\Delta\bm{x}}{2}}\right)}-\varphi{\left({\bm{x^{\prime}}-{\bm{x}_{j}}-\frac{\Delta\bm{x}}{2}}\right)}\right]}\color[rgb]{0,0,0}.$
(28)

The first summation of the final expression (the third last line of the
equation, blue) is for the case of $i=j$, i.e., for the AI group. The second
summation (the last two lines of the equation, brown and red) is for the case
of $i\neq j$, that is, the CI group. Using the Dirac delta function,
${\bm{x}_{p}}$, and ${\bm{x}_{q}}$, the physical spatial coordinates in the
sample, Eq. (28), are rewritten as

$s{\left({\bm{x^{\prime}}+\frac{\Delta\bm{x}}{2}}\right)}s^{*}{\left({\bm{x^{\prime}}-\frac{\Delta\bm{x}}{2}}\right)}\\\
=\color[rgb]{0.2,0.2,0.6}\iint\sum_{i=1}^{N}\delta{\left({{\bm{x}_{p}}-{\bm{x}_{i}}}\right)}B^{2}({\bm{x}_{p}}){{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{p}}+\frac{\Delta\bm{x}}{2}\right)}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{p}}-\frac{\Delta\bm{x}}{2}\right)}\\\
\color[rgb]{0.2,0.2,0.6}\times\exp{i\left[\varphi\left(\bm{x^{\prime}}-{\bm{x}_{p}}+\frac{\Delta\bm{x}}{2}\right)-\varphi\left(\bm{x^{\prime}}-{\bm{x}_{p}}-\frac{\Delta\bm{x}}{2}\right)\right]}d{\bm{x}_{p}}\\\
+\color[rgb]{0.43,0.21,0.1}\iiiint\sum_{i=1}^{N}\sum_{j=1}^{N}\delta{\left({{\bm{x}_{p}}-{\bm{x}_{i}},{\bm{x}_{q}}-{\bm{x}_{j}}}\right)}B({\bm{x}_{p}})B({\bm{x}_{q}})\color[rgb]{0.8,0.0,0.0}e^{i(\phi_{pq})}\\\
\color[rgb]{0.8,0.0,0.0}\times{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{p}}+\frac{\Delta\bm{x}}{2}\right)}{{P_{a}}\left(\bm{x^{\prime}}-{\bm{x}_{q}}-\frac{\Delta\bm{x}}{2}\right)}e^{i\left[\varphi\left(\bm{x^{\prime}}-{\bm{x}_{p}}+\frac{\Delta\bm{x}}{2}\right)-\varphi\left(\bm{x^{\prime}}-{\bm{x}_{q}}-\frac{\Delta\bm{x}}{2}\right)\right]}d{\bm{x}_{p}}d{\bm{x}_{q}}\\\
\color[rgb]{0.43,0.21,0.1}-\mathrm{the\,\,\,first\,\,\,integral\,\,\,of\,\,\,R.H.S.}$
(29)

The subtraction on the last line is because of the same reason as that in Eq.
(7). By substituting the scatterer density maps [Eqs. (8) and (9)] into this
equation, we obtain the continuous form, i.e., Eq. (22).

## References

  * [1] D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, and C. A. Puliafito, “Optical coherence tomography,” Science 254, 1178–1181 (1991).
  * [2] C. A. Curcio, J. D. Messinger, K. R. Sloan, A. Mitra, G. McGwin, and R. F. Spaide, “Human chorioretinal layer thicknesses measured in macula-wide, high-resolution histologic sections,” Invest. Ophthalmol. Vis. Sci. 52, 3943–3954 (2011).
  * [3] E. Sattler, R. Kästle, and J. Welzel, “Optical coherence tomography in dermatology,” J. Biomed. Opt. 18, 061224 (2013).
  * [4] J. M. Schmitt, A. R. Knuettel, A. H. Gandjbakhche, and R. F. Bonner, “Optical characterization of dense tissues using low-coherence interferometry,” in _Holography, Interferometry, and Optical Pattern Recognition in Biomedicine III,_ vol. 1889 H. P. M.D., ed., International Society for Optics and Photonics (SPIE, 1993), pp. 197 – 211.
  * [5] J. M. Schmitt and A. Knüttel, “Model of optical coherence tomography of heterogeneous tissue,” J. Opt. Soc. Am. A 14, 1231–1242 (1997).
  * [6] L. Thrane, H. T. Yura, and P. E. Andersen, “Analysis of optical coherence tomography systems based on the extended Huygens–Fresnel principle,” J. Opt. Soc. Am. A 17, 484–490 (2000).
  * [7] D. J. Smithies, T. Lindmo, Z. Chen, J. S. Nelson, and T. E. Milner, “Signal attenuation and localization in optical coherence tomography studied by Monte Carlo simulation,” Phys. Med. Biol. 43, 3025 (1998).
  * [8] G. Yao and L. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. 44, 2307 (1999).
  * [9] P. R. T. Munro, A. Curatolo, and D. D. Sampson, “Full wave model of image formation in optical coherence tomography applicable to general samples,” Opt. Express 23, 2541–2556 (2015).
  * [10] T. Brenner, D. Reitzle, and A. Kienle, “Optical coherence tomography images simulated with an analytical solution of Maxwell’s equations for cylinder scattering,” J. Biomed. Opt. 21, 45001 (2016).
  * [11] P. R. T. Munro, “Three-dimensional full wave model of image formation in optical coherence tomography,” Opt. Express 24, 27016–27031 (2016).
  * [12] K. C. Zhou, R. Qian, A.-H. Dhalla, S. Farsiu, and J. A. Izatt, “Unified $k$-space theory of optical coherence tomography,” Adv. Opt. Photon. 13, 462–514 (2021).
  * [13] R. D. Allen, G. B. David, and G. Nomarski, “The zeiss-Nomarski differential interference equipment for transmitted-light microscopy,” Z. Wiss. Mikrosk. 69, 193–221 (1969). Zeitschrift fur wissenschaftliche Mikroskopie und mikroskopische Technik.
  * [14] D. Oida, K. Tomita, K. Oikawa, T.-A. Wang, S. Makita, M.-T. Tsai, and Y. Yasuno, “Computational multi-directional optical coherence tomography for visualizing the microstructural directionality of the tissue,” Biomed. Opt. Express 12, 3851–3864 (2021).
  * [15] D. Oida, K. Oikawa, T.-A. Wang, M.-T. Tsai, S. Makita, and Y. Yasuno, “Computational multi-directional optical coherence tomography for visualizing the microstructural directionality of the tissue,” Proc. SPIE 11630, 116301C (2021).
  * [16] M. Born and E. Wolf, _Principles of optics: 60th anniversary edition_ (Cambridge University Press, 2019), chap. 13. Scattering from inhomogeneous media, pp. 695–734, 7th ed.
  * [17] T. S. Ralston, D. L. Marks, F. Kamalabadi, and S. A. Boppart, “Deconvolution methods for mitigation of transverse blurring in optical coherence tomography,” IEEE Trans. Image Process. 14, 1254–1264 (2005).
  * [18] K. W. Gossage, C. M. Smith, E. M. Kanter, L. P. Hariri, A. L. Stone, J. J. Rodriguez, S. K. Williams, and J. K. Barton, “Texture analysis of speckle in optical coherence tomography images of tissue phantoms,” Phys. Med. Biol. 51, 1563 (2006).
  * [19] C. J. Cogswell and C. J. R. Sheppard, “Confocal differential interference contrast (DIC) microscopy: including a theoretical analysis of conventional and confocal DIC imaging,” J. Microsc. 165, 81–101 (1992).
  * [20] J. Canny, “A computational approach to edge detection,” IEEE Trans. Pattern Anal. Mach. Intell. PAMI-8, 679–698 (1986).
  * [21] K. Oikawa, D. Oida, S. Makita, and Y. Yasuno, “Bulk-phase-error correction for phase-sensitive signal processing of optical coherence tomography,” Biomed. Opt. Express 11, 5886–5902 (2020).
  * [22] Y. Yasuno, J. ichiro Sugisaka, Y. Sando, Y. Nakamura, S. Makita, M. Itoh, and T. Yatagai, “Non-iterative numerical method for laterally superresolving Fourier domain optical coherence tomography,” Opt. Express 14, 1006–1020 (2006).
  * [23] L. Zhu, S. Makita, D. Oida, A. Miyazawa, K. Oikawa, P. Mukherjee, A. Lichtenegger, M. Distel, and Y. Yasuno, “Computational refocusing of Jones matrix polarization-sensitive optical coherence tomography and investigation of defocus-induced polarization artifacts,” Biomed. Opt. Express 13, 2975–2994 (2022).
  * [24] M. Villiger and T. Lasser, “Image formation and tomogram reconstruction in optical coherence microscopy,” J. Opt. Soc. Am. A 27, 2216–2228 (2010).
  * [25] T. S. Ralston, D. L. Marks, P. S. Carney, and S. A. Boppart, “Inverse scattering for optical coherence tomography,” J. Opt. Soc. Am. A 23, 1027–1037 (2006).
  * [26] Y. Wang, B. A. Bower, J. A. Izatt, O. Tan, and D. Huang, “In vivo total retinal blood flow measurement by Fourier domain Doppler optical coherence tomography,” J. Biomed. Opt. 12, 041215 (2007).
  * [27] B. R. White, M. C. Pierce, N. Nassif, B. Cense, B. H. Park, G. J. Tearney, B. E. Bouma, T. C. Chen, and J. F. de Boer, “In vivo dynamic human retinal blood flow imaging using ultra-high-speed spectral domain optical Doppler tomography,” Opt. Express 11, 3490–3497 (2003).
  * [28] R. A. Leitgeb, L. Schmetterer, W. Drexler, A. F. Fercher, R. J. Zawadzki, and T. Bajraszewski, “Real-time assessment of retinal blood flow with ultrafast acquisition by color Doppler Fourier domain optical coherence tomography,” Opt. Express 11, 3116–3121 (2003).
  * [29] U. Tricoli and R. Carminati, “Modeling of full-field optical coherence tomography in scattering media,” J. Opt. Soc. Am. A 36, C122–C129 (2019).
  * [30] V. Barolle, J. Scholler, P. Mecê, J.-M. Chassot, K. Groux, M. Fink, A. C. Boccara, and A. Aubry, “Manifestation of aberrations in full-field optical coherence tomography,” Opt. Express 29, 22044–22065 (2021).


# Low rank approximation of positive semi-definite symmetric matrices using
Gaussian elimination and volume sampling

M. Hegland  F. de Hoog

(1 December 2020)

###### Abstract

Positive semi-definite matrices commonly occur as normal matrices of least
squares problems in statistics or as kernel matrices in machine learning and
approximation theory. They are typically large and dense. Thus algorithms to
solve systems with such a matrix can be very costly. A core idea to reduce
computational complexity is to approximate the matrix by one with a low rank.
The optimal and well understood choice is based on the eigenvalue
decomposition of the matrix. Unfortunately, this is computationally very
expensive.

Cheaper methods are based on Gaussian elimination but they require pivoting.
We will show how invariant matrix theory provides explicit error formulas for
an averaged error based on volume sampling. The formula leads to ratios of
elementary symmetric polynomials on the eigenvalues. We discuss some new an
old bounds and include several examples where an expected error norm can be
computed exactly.

##  1  Introduction

The eigenvalue decomposition of a _real symmetric semi-positive definite
matrix $M\in\mathbb{R}^{n,n}$ of rank $r$_ is

$M=Q\Lambda Q^{T}$ (1)

where the factor $Q\in\mathbb{R}^{n,r}$ has $r$ orthonormal columns and
$\Lambda\in\mathbb{R}^{r,r}$ is a diagonal matrix with elements
$\lambda_{1}\geq\cdots\geq\lambda_{r}>0$. Such matrices are common in
applications in machine learning and information retrieval among others [4].
The normal matrix $M=X^{T}X$ which occurs in least squares problems is one
instance.

Often $M$ is very large, dense and unstructured. If the rank $r$ is small,
however, the eigenvalue decomposition (1) shows that it is possible to
represent $M$ using the $nr+r$ matrix elements of the matrices $Q$ and
$\Lambda$. If many of the eigenvalues $\lambda_{i}\geq 0$ are close to zero
one may set these values to zero. So given a matrix $M$ one then might compute
the eigenvalue decomposition, remove small eigenvalues and their corresponding
eigenvectors from the factors $\Lambda$ and $Q$ to get a good approximation
$M_{k}$ which is represented by $nk+k$ real numbers. In fact, the Eckart-
Young-Mirsky theorem [1] states that the approximation $M_{k}$ is optimal.
Such an approximation leads to improved computational performance when used in
algorithms requiring matrix vector products, leads to higher stability and is
useful in data analysis (see principle component analysis). However, the
computational cost of the eigenvalue decomposition is typically of order
$O(n^{3})$ requiring the storage of $n^{2}$ numbers. Thus for large $n$, this
approach is often not feasible in practice.

Thus there is a real need to have a faster algorithm which obtains close to
optimal approximation of $M$. Here we consider a popular example of such an
algorithm which is often referred to CUR or pseudo-skeleton approximation [3].
This algorithm selects $k$ columns (or rows) of the symmetric matrix $M$ and
uses them to approximate $M$. More specifically one may represent $M$ by the
equation

$M=P^{T}\begin{bmatrix}A&B^{T}\\\ B&C\end{bmatrix}P^{T}$ (2)

where $P$ permutes the columns such that the selected ones are moved to the
top and $A\in\mathbb{R}^{k,k}$. The CUR approximation then is of the form

$\widehat{M}_{k}=P^{T}\begin{bmatrix}A&B^{T}\\\ B&BA^{-1}B^{T}\end{bmatrix}P.$
(3)

One can show that this approximation is of rank $k$ and in the following
sections we will investigate the error of this approximation. For the CUR
approximation (3) to be defined, $A$ needs to be invertible. One approach is
to select the $k$ columns for which the determinant of $A$ is maximal. If the
matrix $M$ has a rank $r\geq k$ this choice would guarantee that $A$ is
invertible. Here we consider an approach which selects the $k$ columns at
random with probability proportional to the determinant of $A$. This choice
has been termed _volume sampling_ [2]. In this case the probability of
selecting $k$ columns which lead to a non-invertible $A$ is zero. For this
method one can get an exact expression for the expectation of a suitable norm
of the error. This is similar to the optimal case. In contrast to the optimal
case, however, the determination of the error is computationally often not
feasible even if all the eigenvalues are known. Here we will study this
further and present upper bounds for these expected errors.

One suitable norm for the error analysis is the nuclear norm (also called
trace norm or Schatten 1 norm). The nuclear norm of $M$ is the sum of its
singular values. In the case of semi-definite symmetric matrices $M$ the
singular values are equal to the eigenvalues and one has

$\lVert M\rVert_{*}=\sum_{i=1}^{n}\lambda_{i}$ (4)

where $\lVert M\rVert_{*}$ denotes the nuclear norm. As we will always assume
that the eigenvalues are numbered in decreasing order
($\lambda_{i+1}\leq\lambda_{i}$), the nuclear norm of the error of the optimal
approximation is

$\lVert M_{k}-M\rVert_{*}=\sum_{i=k+1}^{n}\lambda_{i}.$ (5)

In Section 2 the volume sampling CUR approach is discussed and a formula for
the expected error in terms of matrix invariants is established. The errors as
functions of the eigenvalues are further discussed in Section 3 and two
special types of matrices are considered in more depth.

##  2  Rank k approximation and expected error in terms of matrix invariants

In order to establish the framework for volume sampling we introduce the
sample space to be the symmetric group $\Omega=S_{n}$. Then any sample
$\omega\in S_{n}$ is a permutation of a set with $n$ elements. The symmetric
group is the structure of the set of permutation matrices of arrays with $n$
elements and we will denote the permutation matrix defined by some
$\omega\in\Omega$ as $P_{\omega}$. The function mapping $\omega$ to the set of
matrices defined by $M(\omega)=P_{\omega}MP_{\omega}^{T}$ is then a matrix
valued random variable. We will denote the blocks defined in equation (2) of
$M(\omega)$ by $A(\omega)$, $B(\omega)$ and $C(\omega)$ and for simplicity
denote the corresponding variables by $M$, $A$, $B$, and $C$. Finally, we
define the probability of some $\omega\in\Omega$ to be

$p_{k}(\omega):=\frac{\operatorname{det}A(\omega)}{\sum_{\xi\in\Omega}\operatorname{det}A(\xi)}.$
(6)

With this framework we can now define the expected error of the CUR
approximation to be

$\operatorname{\mathcal{E}}\left(\lVert\widehat{M}_{k}-M\rVert_{*}\right):=\sum_{p_{k}(\omega)_{\neq}0}p_{k}(\omega)\,\lVert
C(\omega)-B(\omega)A(\omega)^{-1}B(\omega)^{T}\lVert_{*}$ (7)

i.e. the expectation is a sum over permutations $\omega$ with nonzero
probability.

The determinant $\operatorname{det}A(\gamma)$ is a $k$-th order principle
minor of the unpermuted matrix $M$. On can see that each principle minor of
$M$ occurs $k!(n-k)!$ times when cycling through all the elements of
$\gamma\in\Omega$. Let $c_{j}(M)$ denote the sum all $j$-th principle minors
of $M$ in the following theorem.

######  Theorem 1 (expected error of volume sampling CUR).

$\operatorname{\mathcal{E}}\left(\lVert\widehat{M}_{k}-M\rVert_{*}\right)=(k+1)\frac{c_{k+1}(M)}{c_{k}(M)}$
(8)

In the proof we will use the lemma:

######  Lemma 2.

If $\begin{bmatrix}A&b\\\ b^{T}&\gamma\end{bmatrix}$ is positive semidefinite
and $\operatorname{det}A=0$ then

$\operatorname{det}\begin{bmatrix}A&b\\\ b^{T}&\gamma\end{bmatrix}=0.$

#### Proof:

As $\operatorname{det}A=0$ there exists an $x\neq 0$ such that $Ax=0$. From
the semi positive-definiteness one then gets for all $\eta\in\mathbb{R}$

$\begin{bmatrix}x^{T}&\eta\end{bmatrix}\begin{bmatrix}A&b\\\
b^{T}&\gamma\end{bmatrix}\begin{bmatrix}x\\\ \eta\end{bmatrix}=2\eta
x^{T}b+\eta^{2}\gamma\geq 0$

and thus $b^{T}x=0$. Thus $\begin{bmatrix}x\\\ 0\end{bmatrix}$ is in the null
space of $\begin{bmatrix}A&b\\\ b^{T}&\gamma\end{bmatrix}$.  $\spadesuit$

#### Proof of Theorem 1:

As each minor occurs $k!(n-k)!$ times in the sequence $A(\Omega)$ one has

$\sum_{\omega\in\Omega}detA(\omega)=k!(n-k)!c_{k}(M)$

and consequently
$p_{k}(\omega)=\frac{\operatorname{det}A(\omega)}{k!(n-k)!c_{k}(M)}.$ As
$C-BA^{-1}B^{T}$ is positive semi-definite one gets

$\lVert
BA^{-1}B^{T}-C\rVert_{*}=\operatorname{trace}C-BA^{-1}B^{T}=\sum_{i=1}^{n-k}c_{i,i}-b_{i}^{T}A^{-1}b_{i}=\frac{\sum_{i=1}^{n-k}\operatorname{det}\begin{bmatrix}A&b_{i}\\\
b_{i}^{T}&c_{i,i}\end{bmatrix}}{\operatorname{det}A}$

where $b_{i}^{T}$ is the $i$-th row of $B$. Multiplying both sides with
$\operatorname{det}A$ then summing over all permutations and applying Lemma 2
gives

$\sum_{p(\omega)\neq 0}(\operatorname{det}A(\omega))\lVert
B(\omega)A(\omega)^{-1}B(\omega)^{T}-C(\omega)\lVert_{*}=\sum_{\omega\in\Omega}\sum_{i=1}^{n-k}\
det\begin{bmatrix}A&b_{i}\\\ b_{i}^{T}&c_{i,i}\end{bmatrix}$

Interchanging the order of the double sum gives $(k+1)!(n-k)!c_{k+1}(M)$ and
inserting $p_{k}(A)=(\operatorname{det}A)/(k!(n-k)!c_{k}(M))$ from above
completes the proof.  $\spadesuit$

##  3 The case of diagonal matrices or error and stability bounds in terms of
the eigenvalues

The matrix invariants $c_{k}(M)$ show up as coefficients of the characteristic
polynomial (up to signs). Thus they are also invariant under orthogonal
similarity transforms and thus

$c_{k}(M)=c_{k}(\Lambda)=e_{k}(\lambda_{1},\ldots,\lambda_{n})$

where $e_{k}$ are the _elementary symmetric polynomials_. The evaluation of
elementary symmetric polynomials is computationally very demanding in general.
We will thus focus on upper bounds. A simple bound is given in Proposition 3.
This bound is tight and gives a good indication of what the error is for
spectra which are rapidly decreasing. For spectra like $\lambda_{i}=1/i^{2}$,
however, this bound is highly overestimating the error. In this section we
will discuss ways to get better bounds for these cases.

$0$$250$$500$$750$$1000$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$i$$\lambda_{i}$spectrum
and upper bound$\lambda_{i}=1/i^{2}$piecewise upper bound

$0$$10$$20$$30$$40$$50$$60$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$k$e_{k+1}/e_{k}$upper
bound using pw constant spectrum$\lambda_{i}=1/i^{2}$piecewise upper
boundcurrent bound

Figure 1:  Eigenvalues and the ratios $e_{k+1}/e_{k}$

First we demonstrate the performance of these methods with a computational
experiment. Let $\lambda_{i}=1/i^{2}$ for $i=1,\ldots,n$. We then define a
piecewise constant array $\mu$ by

$\mu_{i}=\frac{1}{4^{l}},\quad\text{where $i=2^{l}+j$, $j=1,\ldots,2^{l}-1$
and $l=0,\ldots,l_{\text{max}}-1$}.$

One can show that $\mu_{i}\geq\lambda_{i}$ for $i=1,\ldots,n$. From the
monotonicity (Lemma 3) one sees that

$\frac{e_{k+1}(\lambda)}{e_{k}(\lambda)}\leq\frac{e_{k+1}(\mu)}{e_{k}(\mu)}.$

This is illustrated computationally in Figure 1 where the simple bound from
Proposition 3 is also displayed to demonstrate the superiority of the new
bound. While not as simple as the simple bound, the bound we give here is
substantially more accurate and useful for error bounds. In contrast to the
values for $\lambda$ the upper bound can be computed even for very large $n$
and medium sized $k$ as explained below.

An upper bound for $e_{k+1}/e_{k}$ which does not require any properties of
the eigenvalues is obtained from

######  Proposition 3 (simple bound).

Let $n,k\in\mathbb{N}$, $n>k$ and $\lambda_{i}\geq 0$ for all $i=1,\ldots,n$.
Then

$e_{k+1}(\lambda_{1},\ldots,\lambda_{n})\leq
e_{k}(\lambda_{1},\ldots,\lambda_{n})\,e_{1}(\lambda_{k+1},\ldots,\lambda_{n}).$

#### Proof:

The proof is by induction where the main step is (the second line uses the
induction assumption)

$\displaystyle e_{k+1}(\lambda_{1},\ldots,\lambda_{n+1})$
$\displaystyle=e_{k}(\lambda_{1},\ldots,\lambda_{n})\,\lambda_{n+1}+e_{k+1}(\lambda_{1},\ldots,\lambda_{n})$
$\displaystyle\leq
e_{k}(\lambda_{1},\ldots,\lambda_{n})\,\lambda_{n+1}+e_{k}(\lambda_{1},\ldots,\lambda_{n})\,e_{1}(\lambda_{k+1},\ldots,\lambda_{n})$
$\displaystyle=e_{k}(\lambda_{1},\ldots,\lambda_{n})\,e_{1}(\lambda_{k+1},\ldots,\lambda_{n+1}).$

$\spadesuit$ Thus $e_{k+1}/e_{k}$ is bounded by the optimal error and by
Theorem 1 the CUR error is bounded by $k+1$ times the optimal error.

In the case where $\lambda_{i}=q^{i}$ for some $q\in(0,1)$ one can compute the
elementary symmetric polynomials explicitly as

######  Proposition 4 (power eigenvalues).

$e_{k}(1,q,\ldots,q^{n-1})=q^{k(k-1)/2}\,\prod_{i=1}^{k}\frac{1-q^{n-i+1}}{1-q^{i}},\quad\text{for
all $k,n\in\mathbb{N}$}.$ (9)

#### Proof:

Let $e_{k,n}:=e_{k}(1,q,\ldots,q^{n-1})$. We then use induction over $n$ to
show that equation (9) holds for all $k$.

First, for $n=1$ one has $e_{0,1}=e_{1,1}=1$ and $e_{k,1}=0$ by definition of
$e_{k}$. Thus the claimed result holds for $n=1$.

We now have to show that if equation (9) holds for some $n$ and all $k$ it
also holds when $n$ is replaced by $n+1$. For this we use the following
recursion for the symmetric elementary polynomials:

$e_{k,n+1}=q^{n}e_{k-1,n}+e_{k,n}.$

Then one verifies that the claimed equation (9) is equivalent to

$e_{k,n}=\prod_{i=1}^{k}\frac{q^{i-1}-q^{n}}{1-q^{i}}.$

It follows that

$e_{k,n}=\frac{q^{k-1}-q^{n}}{1-q^{k}}\,e_{k-1,n}.$

One then gets

$\displaystyle e_{k,n+1}$
$\displaystyle=q^{n}e_{k-1,n}+e_{k,n}=\left(q^{n}+\frac{q^{k-1}-q^{n}}{1-q^{k}}\right)\,e_{k-1,n}$
$\displaystyle=\frac{q^{k-1}(1-q^{n+1})}{1-q^{k}}\,e_{k-1,n}=\frac{1-q^{n+1}}{1-q^{k}}\,\prod_{i=1}^{k-1}\frac{q^{i}-q^{n+1}}{1-q^{i}}$
$\displaystyle=\prod_{i=1}^{k}\frac{q^{i-1}-q^{n+1}}{1-q^{i}}.$

$\spadesuit$ The ratio $e_{k+1}/e_{k}=\frac{q^{k}-q^{n}}{1-q^{k+1}}$ obtained
from this result then leads to the expected error

$(k+1)\frac{e_{k+1}}{e_{k}}=\frac{k+1}{1+q+\cdots+q^{k}}(q^{k}+\cdots+q^{n}).$

One sees that the CUR method is competitive for small $q$ but not for
$q\approx 1$.

We now derive some results which are useful in the derivation and computation
of error bounds. First we pad the eigenvalue vectors with zeros so that
$\lambda\in\ell_{0}$, the set of series which are nonzero for finitely many
indices.

For any $k$ we define the head of $\lambda$ to be
$\lambda^{h}=(\lambda_{1},\ldots,\lambda_{k},0,\ldots)$ and the tail
$\lambda^{t}=(\lambda_{k+1},\ldots,\lambda_{n},0,\ldots)$. We introduce a
concatenation of two spectra $\lambda$ and $\mu$ denoted by $[\lambda,\mu]$
where

$[\lambda,\mu]=(\lambda_{1},\ldots,\lambda_{n},\mu_{1},\ldots,\mu_{m},0,\ldots)$

which potentially is reordered for size but in the cases considered here we
have $\mu_{1}\leq\lambda_{n}$. Thus $\lambda=[\lambda^{h},\lambda^{t}]$. The
decrease of the tail is modelled by the sequence $\rho$ with

$\rho_{i}=\lambda^{t}/\lambda_{k+1}.$

We then introduce a function
$f:\mathbb{R}_{0}^{\infty}\rightarrow\mathbb{R}^{k+2}$ with components

$f_{i}(\lambda)=e_{i}(\lambda).$

We now introduce the _convolution_ of two elements of $\mathbb{R}_{+}^{m}$ by

$(u*v)_{i}=\sum_{i+j=k}u_{i}v_{j}.$

As the $e_{i}$ are coefficients of a characteristic polynomial one has

######  Lemma 5 (convolution theorem).

$f([\lambda,\mu])=f(\lambda)*f(\mu).$

The next lemma is a consequence of the fact that $f_{i}$ is a homogeneous
$i-1$st degree polynomial.

######  Lemma 6 (scaling lemma).

The $i$-th component $f_{i}$ of $f$ satisfies

$f_{i}(s\lambda)=s^{i-1}f_{i}(\lambda).$

One then has

######  Proposition 7.

Let $u^{h}=f(\lambda^{h})$ and $w=f(\rho)$. Then

$\frac{e_{k+1}(\lambda)}{e_{k}(\lambda)}=\gamma\lambda_{k+1},$

where

$\gamma=\frac{\sum_{i=0}^{k}\lambda_{k+1}^{k-i}\;w_{k+1-i}u^{h}_{i}}{\sum_{i=0}^{k}\;\lambda_{k+1}^{k-i}w_{k-i}u^{h}_{i}}.$

#### Proof:

By the convolution theorem and the definition of $\rho$ one has

$u=u^{t}*u^{h}$

where $u=f(\lambda)$ and $u^{t}=f(\lambda^{t})=f(\lambda_{k+1}\rho)$.

Using the scaling lemma and $u^{h}_{k+1}=0$ one then gets

$\frac{u_{k+1}}{u_{k}}=\lambda_{k+1}\frac{\sum_{i=0}^{k}\lambda_{k+1}^{k-i}\;w_{k+1-i}u^{h}_{i}}{\sum_{i=0}^{k}\;\lambda_{k+1}^{k-i}w_{k-i}u^{h}_{i}}.$

$\spadesuit$

We will now show results used to obtain bounds for the case of slowly
decreasing $\rho_{i}$. In these cases one observes that the sequence
$w_{i+1}/w_{i}$ first increases before it decreases. The main tool to obtain
bounds is the monotonicity in of the ratios $e_{k+1}(\lambda)/e_{k}(\lambda)$
in $\lambda$.

######  Lemma 8 (monotonicity).

Let $0<\lambda\leq\mu$ component wise then

$\frac{e_{k+1}(\lambda)}{e_{k}(\lambda)}\leq\frac{e_{k+1}(\mu)}{e_{k}(\mu)}.$

This lemma can be proven directly but is also a consequence of a result by
Marcus and Lopes [5]:

$\frac{e_{k+1}(\lambda+\mu)}{e_{k}(\lambda+\mu)}\geq\frac{e_{k+1}(\lambda)}{e_{k}(\lambda)}+\frac{e_{k+1}(\mu)}{e_{k}(\mu)}.$

This holds for $\mu$ and $\lambda$ being nonegative.

A direct application of the convolution theorem gives the representation

######  Proposition 9.

$f(\lambda)={u^{l_{\text{max}}}}*\cdots*u^{0}$ (10)

where $u^{l}=f(q^{l}\epsilon_{2^{l}})$ has the components

$u^{l}_{j}=q^{l(j-1)}\binom{2^{k}}{j-1},\quad j=1,\ldots,2^{l}.$

This proposition is used to show that the ratio $e_{k+1}/e_{k}$ is
computationally feasible for the piecewise spectrum $\mu$ used in Figure 1.
The complexity of computing the $l_{\text{max}}$ convolutions of size $k+1$ is
of order $O(l_{\text{max}}(k+1)^{2})$. Using this formula is typically much
faster than using the standard recursions which results in a total complexity
of $O(n(k+1))$ for very large $n$. One observes that often only a small number
of components of the $u^{j}$ are substantially different from zero so that the
complexity can be further reduced.

##  4 Conclusion

While the CUR method combined with volume sampling admits an explicit and
exact formula for the expected approximation error of the computed rank $k$
approximation there is little known about the performance of the method for
slowly decreasing eigenvalues. This paper provides new error bounds for the
CUR method with volume sampling which shows that this approach is competitive
with the optimal approach for slowly decreasing eigenvalues even for large
approximation ranks.

reviews the basic theory and provides a framework and some simple examples
which provide some insights on how the CUR method performs well in particular
for these cases of slowly decreasing eigenvalues. Future work may consider the
effect of the initial eigenvalues and the rate of decrease of the tail of the
spectrum in more detail using the convolution formula provided here.

## References

  * [1] Achiya Dax  “On extremum properties of orthogonal quotients matrices”  In _Linear Algebra and its Applications_ 432.5, 2010, pp. 1234–1257  DOI: 10.1016/j.laa.2009.10.034
  * [2] Amit Deshpande, Luis Rademacher, Santosh Vempala and Grant Wang  “Matrix Approximation and Projective Clustering via Volume Sampling”  In _Proceedings of the Seventeenth Annual ACM-SIAM Symposium on Discrete Algorithm_ , SODA ’06  Miami, Florida: Society for IndustrialApplied Mathematics, 2006, pp. 1117–1126
  * [3] S.A. Goreinov, E.E. Tyrtyshnikov and N.L. Zamarashkin  “A theory of pseudoskeleton approximations”  In _Linear Algebra and its Applications_ 261.1, 1997, pp. 1 –21  DOI: https://doi.org/10.1016/S0024-3795(96)00301-1
  * [4] Michael W. Mahoney and Petros Drineas  “CUR matrix decompositions for improved data analysis”  In _Proceedings of the National Academy of Sciences_ 106.3, 2009, pp. 697–702  DOI: 10.1073/pnas.0803205106
  * [5] M. Marcus and L. Lopes  “Inequalities for Symmetric Functions and Hermitian Matrices”  In _Canadian Journal of Mathematics_ 9 Cambridge University Press, 1957, pp. 305–312  DOI: 10.4153/CJM-1957-037-9

### Author address