# Application of Neural Network Machine Learning to Solution of Black-Scholes
Equations

Mikhail V. Klibanov, Kirill V. Golubnichiy and Andrey V. Nikitin  Department
of Mathematics and Statistics  Department of Mathematics  University of North
Carolina at Charlotte, Charlotte, USA  University of Washington, Seattle, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper presents a novel way to predict options price for one day in
advance, utilizing the method of Quasi-Reversibility for solving the Black-
Scholes equation. The Black-Scholes equation is solved forwards in time, which
is an ill-posed problem. Thus, Tikhonov regularization via the Quasi-
Reversibility Method is applied. This procedure allows to forecast stock
option prices for one trading day ahead of the current one. To enhance these
results, the Neural Network Machine Learning is applied on the second stage.
Real market data are used. Results of Quasi-Reversibility Method and Machine
Learning method are compared in terms of accuracy, precision and recall.

Keywords:

The Black-Scholes equation, Ill-posed problem, regularization method,
parabolic equation with the reversed time, Machine Learning, neural network.

##  1 Introduction

This paper discusses a new empirical mathematical model for generating more
accurate option trading strategy using initial and boundary conditions for the
underlying stock. The idea was initially proposed in [5]. The basis for this
idea is the Black-Scholes equation. In mathematical finance, the Black-Scholes
equation is a parabolic partial differential equation that determines the
dynamics of the price of European options [7].

The time at a given time $t$ will occur is $\tau,$

$\tau=T-t.$ (1.1)

$f(s)$ be the payoff function of that option at the maturity time $t=T$ and
$s$ is the stock price. Let’s assume that the risk-free interest rate equals
zero. The function $u\left(s,\tau\right)$ is the price of that option and the
variable $\tau$ is the one defined in (1.1). Let’s assume that this function
$u\left(s,\tau\right)$ satisfies the Black-Scholes equation with the
volatility coefficient ${\sigma}$ [11, Chapter 7, Theorem 7.7]:

$\begin{split}&\frac{\partial
u(s,\tau)}{\partial\tau}=\frac{{\sigma}^{2}}{2}s^{2}\frac{\partial^{2}u(s,\tau)}{\partial
s^{2}},\\\ &u(s,0)=f(s),\end{split}$ (1.2)

The payoff function is $f(s)=\max\left(s-K,0\right)$, where $K$ is the strike
price [11] and $s>0.$

The option price function is defined by the Black-Scholes formula:

$u(s.\tau)=s\Phi(\theta_{+}(s,\tau))-e^{-r\tau}K\Phi(\theta_{-}(s,\tau)),$
(1.3)

Based on the Itô formula, we have:

$du=(-\frac{\partial
u(s,T-t)}{\partial\tau}+\frac{{\sigma}^{2}}{2}s^{2}\frac{\partial^{2}u(s,T-t)}{\partial
s^{2}})dt+\sigma s\frac{\partial u(s,T-t)}{\partial s}dW.$ (1.4)

If equation (1.2) is solved forwards in time to forecast prices of stock
options is an ill-posed inverse problem. For this reason, we used the Method
of Quasi- Reversibility (QRM) that is a version of the Tikhonov regularization
method. Uniqueness, stability and convergence theorems for this method were
formulated in [5] and [3], also, see [10] for proofs.

We have four major questions that we raise in this paper:

  1. 1.

What is the forecast interval of the options prices?

  2. 2.

What are the boundary and initial conditions on the interval for the Black-
Scholes equation?

  3. 3.

What are the values of the volatility coefficient in the future?

  4. 4.

How to solve the Black-Scholes equation forwards in time $t$?

The first three questions are addressed in our new mathematical model. We use
the regularization method of [4] to address the fourth question. Theorems
about stability and convergence of this method are formulated. These theorems
were proven in [4] for a general parabolic equation of the second order where
the main key of this method is based on the method of Carleman estimates.

Let the function $f$ $\in L_{2}(0,\pi)$ and let $T=const.>0.$ To demonstrate
that our problem is ill-posed, we consider the example of based on the problem
for the heat equation with the reversed time

$u_{t}+u_{xx}=0,\left(x,t\right)\in\left(0,\pi\right)\times\left(0,T\right),$
(1.5)

with Dirichlet boundary conditions

$u\left(0,t\right)=0,u\left(\pi,t\right)=0,$ (1.6)

and the initial condition

$u\left(x,0\right)=f\left(x\right).$ (1.7)

The unique solution of this problem is:

$u(x,t)=\sum_{n=1}^{\infty}f_{n}\sin{(nx)}e^{n^{2}t}.$

Consider

$u_{N}(x,t)=\sum_{n=1}^{N}f_{n}\sin{(nx)}e^{n^{2}t}.$

Then

$||u_{N}(x,T)||_{L_{2}(0,\pi)}^{2}=\sum_{n=1}^{N}f_{n}^{2}e^{2n^{2}T}\approx
f_{N}^{2}e^{2N^{2}T}\rightarrow\infty$

as $N,T\rightarrow\infty.$ Hence, the problem (1.5)-(1.7) is severely
unstable.

We conclude therefore that to obtain a more or less accurate solution of the
Black-Scholes equation forwards in time, we need to solve it on a short time
interval $\left(0,T\right)$. To get better accuracy, the regularization method
works only for a short time interval.

Section 3 presents our mathematical model, the method of Quasi Reversibility
as well as the trading strategy. The Quasi Reversibility Method is based on
the minimization of a Tikhonov-like functional $J_{\beta}(u).$ We do this
using conjugate gradient method. The minimization process was performed by
Hyak Next Generation Supercomputer of the research computing club of
University of Washington. The code was parallelized in order to maximize the
performance on supercomputer clusters.

The historical data for stock options was collected from the Bloomberg
terminal [2] of University of Washington. From this data, we obtained about
177,000 minimizers.

Due to ill-posedness of the problem the solution is very sensitive to the
noise in the initial data (stock and option prices for the three days
preceding the day of forecast). Given results of the Quasi Reversibility
Method, we apply on the second stage Machine Learning to reduce the
probability of non-profitable trades caused by wrong option price prognosis
because of the noise in input data.

Section 4 is dedicated to application of binary classification and regression
Neural Network Machine Learning.

Sections 5 and 6 present our results and the summary.

Python with the SciPy and Torch modules were used for implementation of the
method of Quasi-Reversibility and Neural Network Machine Learning (binary
classification and regression).

##  2 The new mathematical model and the method of Quasi-Reversibility

Let’s denote $s$ as the stock price, $t$ as the time and $\sigma(t)$ as the
volatility of the option. The historical implied volatility listed on the
market data of [2] is used in our particular case. We assume that
$\sigma=\sigma(t)$ to avoid other historical data for the volatility. Let’s
call $u_{b}(t)$ and $u_{a}(t)$ the bid and ask prices of the options at the
moment of time $t$ and $s_{b}(t)$ and $s_{a}(t)$ the bid and ask prices of the
stock at the moment of time $t.$ It is also known that

$u_{b}(t)<u_{a}(t)$ (2.1)

and

$s_{b}(t)<s_{a}(t)$ (2.2)

Let’s introduce

$f_{s}(t)=\frac{s_{a}(t)}{s_{b}(t)}-1$ (2.3)

and

$f_{u}(t)=\frac{u_{a}(t)}{u_{b}(t)}-1$ (2.4)

Based on real market data we have observed that usually

$0\leq f_{s}(t)\leq 0.003$ (2.5)

and

$0\leq f_{u}(t)\leq 0.27$ (2.6)

The idea is to approximate the Black-Scholes equation solutions

$Lu=u_{t}+\frac{\sigma^{2}\left(t\right)}{2}s^{2}u_{ss}=0,\left(s,t\right)\in\left(s_{b}\left(0\right),s_{a}\left(0\right)\right)\times\left(0,2\tau\right)=X_{2\tau},$
(2.7)

with Dirichlet boundary conditions

$u\left(s_{b},t\right)=u_{b}\left(t\right),u\left(s_{a},t\right)=u_{a}\left(t\right),\text{
}t\in\left[0,2\tau\right],$ (2.8)

and the initial condition

$u\left(s,0\right)=f\left(s\right),\text{ \ \
}s\in\left[s_{b}\left(0\right),s_{a}\left(0\right)\right].$ (2.9)

Where $L$ is the partial differential operator of the Black-Scholes equation.
Based on Bloomerg terminal we used with End of Day Underlying Price Last, End
of Day Underlying Price Bid, End of Day Underlying Price Ask, $t$ is time,
$\sigma\left(t\right)$ is the volatility of the stock option. It was used
Implied Volatility Using Last Trade Price (IVOL).

$u\left(s,t\right)$ is the price of the stock option. End of Day Option Price
Last, End of Day Option Price Bid and End of Day Option Price Ask are the
notation that we applied in our algorithm.

Problem 1. _Find the function_ $u\in H^{2}\left(X_{2\tau}\right)$_satisfying
conditions ( 2.7)-(2.9)._

This problem considers as ill-posed since we solve equation (2.7) forwards in
time.

Remarks 3.1: _We increase here the required smoothness of the solution from_
$H^{2,1}\left(X_{2\tau}\right)$_to_ $H^{2}\left(X_{2\tau}\right).$

Our algorithm based on solving the inverse problem for the Black-Scholes with
reversed time equation has five steps:

Step 1 (Dimensionless variables).

We require to make our equation dimensionless. $s_{b}<s_{a}.$ Let’s denote
$s_{b}=s_{b}(0),$ $s_{a}=s_{a}(0).$ Dimensionless variables were applied
$x,t^{\prime}$ such that

$x=\frac{s-s_{b}}{s_{a}-s_{b}}$ (2.10)

$t^{\prime}=\frac{t}{255}$ (2.11)

and now we can say that $s$ is $x$ and $t$ is $t^{\prime}.$

According to these substitutions, the equation becomes

$Ru=u_{t}+\sigma^{2}\left(t\right)A(x)u_{xx},$ (2.12)

where

$A(x)=\frac{255}{2}\frac{[x(s_{a}-s_{b})+s_{b})]^{2}}{(s_{a}-s_{b})^{2}}$
(2.13)

$X_{2\tau}=\left\\{\left(x,t\right)\in\left(0,1\right)\times\left(0,2\tau\right)\right\\}.$
(2.14)

$u(x,0)=g(x),x\in(0,1)$ (2.15)

$u(0,t)=u_{b}(t),u(1,t)=u_{a}(t).$ (2.16)

And the operator $L$ in (2.7) is the operator $R$

Step 2 (Interpolation and extrapolation).

Our goal is to forecast option price from ’today’ to ’tomorrow’ and ’the day
after tomorrow’. We do have 255 trading days annually. For this reason, let’s
introduce $\tau>0$ as our unit of time for which we want to make our
prediction the option price. Because we predict option prices having the
information of these prices, as well as of other parameters for ’today’,
’yesterday’ and ’the day before yesterday’, we consider $\tau$ is one trading
day. ’One day’ $\tau=1/255.$ ’Today’ $t=0.$ ’Tomorrow’ $t=\tau.$ ’The day
after tomorrow’ $t=2\tau.$ The variable $s-$ is for interval, i.e
$s\in\left[s_{b}\left(0\right),s_{a}\left(0\right)\right].$ We applied the
idea associated with interpolation discrete values of functions $u_{b}(t),$
$u_{a}(t),$ and $\sigma(t)$ between these three points (the day before
yesterday, yesterday and today) and then extrapolation functions $u_{b}(t),$
$u_{a}(t)$ between three points (today, tomorrow and the day after tomorrow).
Where $t=-2\tau$ is ”the day before yesterday”, $t=-\tau$ is ”yesterday” and
$t=0$ is ”today”. We used quadratic polynomials for both approximation and
extrapolation of values of functions. Thus, these three functions $u_{b}(t),$
$u_{a}(t),$ and $\sigma(t)$ was obtained for a small future time interval, i.e
$(0,2\tau).$ ([5]). Where $u_{b}(t),$ $u_{a}(t)$ were applied for boundary
conditions and $\sigma(t)$ is coefficient function for our problem. The
initial condition was set as $u(x,0)=g(x)=x(u_{a}(0)-u_{b}(0))x+u_{b}(0).$
This function is the result of approximation by linear function due to the
fact that the interval between bid and ask prices is relatively small. The
domain was $X_{2\tau}=\\{(x,t):x\in(0,1),t\in(0,2\tau)\\}.$

Step 3 (Statement of the Problem).

Problem 2. _Assume that functions_

$u_{b}\left(t\right),u_{a}\left(t\right)\in
H^{2}\left[0,2\tau\right],\sigma\left(t\right)\in C^{1}\left[0,2\tau\right].$
(2.17)

_Find the solution_ $u\in H^{2}\left(X_{2\tau}\right)$_of the following
initial boundary value problem:_

$Ru=0\text{ in }X_{2\tau},$ (2.18)
$u\left(0,t\right)=u_{b}\left(t\right),u\left(1,t\right)=u_{a}\left(t\right),t\in\left(0,2\tau\right),$
(2.19) $u\left(x,0\right)=g\left(x\right),x\in\left(0,1\right),$ (2.20)

_where the partial differential operator_ $R$_is defined in ( 2.12), the
function _$A\left(x\right)$_is defined in ( 2.13), the initial condition
_$g\left(x\right)$_is defined in ( 2.15), and the domain _$X_{2\tau}$_is
defined in ( 2.14)._

######  Theorem 2.1.

The following problem (2.12)-(2.15) has one solution $u\in
H^{2,1}(X_{2\tau}).$

The proof of this theorem is [3].

Step 4 (Numerical method of solving the problem. Regularization).

Due to the ill-posedness of the problem, we can not say about existence of the
solution. Thus, it was applied the regularization method:

Let’s consider function $F(x,t)=x(u_{a}(t)-u_{b}(t))+u_{b}(t),$ $(x,t)\in
X_{2\tau}.$ This function $F\in H^{2}(X_{2\tau}).$ It follows from (2.15) and
(2.16) that

$F\left(x,0\right)=g\left(x\right),$ (2.21)

$F(0,t)=u_{b}\left(t\right),F(1,t)=u_{a}\left(t\right).$ (2.22)

We used an unbounded differential operator $R:H^{2,1}(X_{2\tau})\rightarrow
L^{2}(X_{2\tau}),$ where $H^{2,1}(X_{2\tau})$ is a dense linear set in the
space $L^{2}(X_{2\tau}).$ Where

$Ru=u_{t}+\sigma^{2}\left(t\right)A(x)u_{xx}$ (2.23)

Let’s introduce Tikhonov-like functional as:

$J_{\beta}\left(u\right)=\int_{X_{2\tau}}\left(Ru\right)^{2}dsdt+\beta\left\|u\right\|_{H^{2}\left(X_{2\tau}\right)}^{2},$
(2.24)

where $\beta\in\left(0,1\right)$ is the parameter of regularization. To solve
the problem, we minimized the functional $J_{\beta}\left(u\right)$ on the set

$V=\left\\{u\in
H^{2}\left(X_{2\tau}\right):u\left(0,t\right)=u_{b}\left(t\right),u\left(1,t\right)=u_{a}\left(t\right),u\left(x,0\right)=g\left(x\right)\right\\}.$
(2.25)

Step 5 (Minimization Problem).

Minimization Problem 1.
$J_{\beta}:H^{2}\left(X_{2\tau}\right)\rightarrow\mathbb{R}$_is the
regularization Tikhonov functional._

We have used the converting of our partial derivatives from (2.24) into finite
differences. A finite difference grid was applied to cover the domain
$X_{2\tau}.$ The minimization process was to differentiate our functional
$J_{\beta}(u)$ with respect to the values of the function $u(x,t)$ at each
grid points via conjugate gradient method. The point $u=0$ was used for the
starting point. Based on computational study with simulated data we have
realized that the optimal value of the regularization parameter would be
$\beta=0.01.$

Minimization Problem 1 is a QRM for Problem 2. This is an version of the QRM
for problem (2.18)-(2.20). In section 4 we discuss the theory of this specific
version of the QRM. In particular, Theorem 4.2 of section 4 presents
uniqueness of the solution $u\in H^{2}\left(X_{2\tau}\right)$ of Problem 2 and
implies an estimate of the stability of this solution with respect to the
noise in the data. Theorem 4.3 of section 4 shows existence and uniqueness of
the minimizer $u_{\beta}\in H^{2,1}\left(X_{2\tau}\right)$ of the functional
$J_{\beta}\left(u\right)$ on the set $V$ defined in (2.25). We call such a
minimizer “regularized solution” [9]. Theorem 4.4 estimates convergence rate
of regularized solutions to the exact solution of Problem 2 with the noiseless
data. Such estimates depend on the noise level in the data. All proof of these
theorems are presented in [10].

##  3 Analysis

This section is devoted to convergence analysis for Problem 2 of subsection
3.2. This problem is the initial boundary value problem for parabolic equation
(2.18) with the reversed time. The QRM and convergence analysis for this
problem for a more general parabolic operator in $\mathbb{R}^{n}$ with
arbitrary variable coefficients was proposed in [4]. Then theorems were
presented in [5]. However a stability estimate was not a part of [5], such an
estimate was proven in [4]. The same is true for the convergence theorems of
QRM in [4, 5]. The smallness assumption was lifted in [3] via a new Carleman
estimate. Results of [3] for a 1-D case were significantly modified in this
section. Our computations below on a small time interval
$\left(0,2\tau\right)=\left(0,0.00784\right)$ (see [4, 5], [8, Theorem 1 of
section 2 in Chapter 4] might result in the requirement of even a smaller
length of that interval.

###  3.1 Problem statement

Let’s consider a number $T>0$ and introduce $Q_{T}$ as:

$Q_{T}=\left\\{\left(x,t\right)\in\left(0,1\right)\times\left(0,T\right)\right\\}.$

Consider two numbers $b_{0},b_{1}>0$ and $b_{0}<b_{1}.$ Let the function
$b\left(x,t\right)\in C^{1}\left(\overline{Q}_{T}\right)$ satisfies:

$\text{ }\left\|b\right\|_{C^{1}\left(\overline{Q}_{T}\right)}\leq
b_{1},\text{ }b\left(x,t\right)\geq b_{0}\text{ in }Q_{T}.$ (3.1)

We also have functions $\psi_{0}\left(t\right),\psi_{1}\left(t\right)\in
H^{2}\left(0,T\right).$ In the above case of subsection 3.2,

$T=2\tau,b\left(x,t\right)=\sigma^{2}\left(t\right)A(x),\psi_{0}\left(t\right)=u_{b}\left(t\right),\psi_{1}\left(t\right)=u_{a}\left(t\right).$

We now formulate Problem 3, which is a slight generalization of Problem 2.

Problem 3. _Find a solution_ $v\in H^{2}\left(Q_{T}\right)$_of the following
(IBVP):_

$Nv=v_{t}+b\left(x,t\right)v_{xx}=0\text{ in }Q_{T},$ (3.2)
$v\left(0,t\right)=\psi_{0}\left(t\right),v\left(1,t\right)=\psi_{1}\left(t\right),\text{
}t\in\left(0,T\right),$ (3.3)
$v\left(x,0\right)=z\left(x\right)=\psi_{0}\left(0\right)\left(1-x\right)+\psi_{1}\left(0\right)x,\text{
}x\in\left(0,1\right).$ (3.4)

Remark 4.1. _Because Problem 2 is less general than Problem 3, then this
analysis of converegence for Problem 3 also works for Problem 2._

We use the linear function for $v\left(x,0\right)$ in (3.4) is to simplify the
initial condition in (2.20). Now problem 3 is an IBVP for the parabolic
equation (3.2) with the reversed time. For this reason, the problem can be
considered as ill-posed. Assume that the boundary with a noise of the level
$\nu>0$ in (3.3) are in place. Here $\nu$ is a sufficiently small number, i.e.

$\left\|\psi_{0}-\psi_{0}^{\ast}\right\|_{H^{1}\left(0,T\right)}<\nu,\left\|\psi_{1}-\psi_{1}^{\ast}\right\|_{H^{1}\left(0,T\right)}<\nu,$
(3.5)

where functions $\psi_{0}^{\ast},\psi_{1}^{\ast}\in H^{2}\left(0,T\right)$ are
“ideal” noiseless data. we assume that there exists an exact solution
$v^{\ast}\in H^{2}\left(Q_{T}\right)$ of problem (3.2)-(3.4) with these
noiseless data (based on on the theory of Ill-Posed problems). Below we
present estimates how this noise affects the accuracy of the solution of
Problem 3 and also discuss the convergence rate of numerical solutions
obtained by QRM to the exact one as $\nu\rightarrow 0.$

Let’s introduce the version of functional (2.24):

$I_{\beta}\left(v\right)=\int_{Q_{T}}\left(Nv\right)^{2}dxdt+\beta\left\|v\right\|_{H^{2}\left(Q_{T}\right)}^{2}.$
(3.6)

We also have the set $W\subset H^{2}\left(Q_{T}\right),$

$W=\left\\{v\in
H^{2}\left(Q_{T}\right):v\left(0,t\right)=\psi_{0}\left(t\right),v\left(1,t\right)=\psi_{1}\left(t\right),v\left(x,0\right)=z\left(x\right)\right\\}.$
(3.7)

The solution of Problem 3 is approximate solution by solving the following
problem:

Minimization Problem 2. _Minimize the functional_ $I_{\beta}\left(v\right)$_on
the set_ $W$_given in ( 3.7)._

Minimization Problem 2 is QRM for Problem 3.

###  3.2 Theorems

This subsection presents four theorems for Problem 3. All proofs might be
found in [10]. First, let’s introduce the Carleman Weight Function
$\phi_{\alpha}\left(t\right)$ with $\alpha>2$ for the operator
$\partial_{t}+b\left(x,t\right)\partial_{x}^{2}$ as:

$\phi_{\alpha}\left(t\right)=e^{\left(T+1-t\right)^{\alpha}},\text{
}t\in\left(0,T\right).$ (3.8)

As a result, the function $\phi_{\alpha}\left(t\right)$ is decreasing on
$\left[0,T\right]$, $\phi_{\alpha}^{\prime}\left(t\right)<0,$

$\max_{\left[0,T\right]}\phi_{\alpha}\left(t\right)=\psi_{\alpha}\left(0\right)=e^{\left(T+1\right)^{\alpha}},\text{
}\min_{\left[0,T\right]}\phi_{\alpha}\left(t\right)=\phi_{\alpha}\left(T\right)=e.$
(3.9)

Denote

$H_{0}^{2}\left(Q_{T}\right)=\left\\{u\in
H^{2}\left(Q_{T}\right):u\left(0,t\right)=u\left(1,t\right)=0\right\\}.$
(3.10) $H_{0,0}^{2}\left(Q_{T}\right)=\left\\{u\in
H_{0}^{2}\left(Q_{T}\right):u\left(x,0\right)=0\right\\}.$ (3.11)

Theorem 4.1 (Carleman estimate). _Let the coefficient_ $b\left(x,t\right)$_of
the operator_ $N$_satisfies conditions ( 3.1). Then there exist a sufficiently
large number _$\alpha_{0}=\alpha_{0}\left(T,b_{0},b_{1}\right)>2$_and a
constant_ $C=C\left(T,b_{0},b_{1}\right)>0,$_both depending only on listed
parameters, such that the following Carleman estimate holds for the operator_
$N:$__

$\int_{Q_{T}}\left(Nu\right)^{2}\phi_{\alpha}^{2}dxdt\geq
C\sqrt{\alpha}\int_{Q_{T}}u_{x}^{2}\psi_{\alpha}^{2}dxdt+C\alpha^{2}\int_{Q_{T}}u^{2}\phi_{\alpha}^{2}dxdt$
$-C\sqrt{\alpha}\left\|u\right\|_{H^{2}\left(Q_{T}\right)}^{2}-C\lambda\left(T+1\right)^{\alpha}e^{2\left(T+1\right)^{\alpha}}\left\|u\left(x,0\right)\right\|_{L_{2}\left(0,1\right)}^{2},$
(3.12) $\forall\alpha\geq\alpha_{0},\forall u\in H_{0}^{2}\left(Q_{T}\right).$

Carleman estimate (3.12) is the MAIN TOOL to proofs of Theorems 4.2, 4.4.

Theorem 4.2 (Hölder stability estimate for Problem 3 and uniqueness). _Let the
coefficient_ $b\left(x,t\right)$_of the operator_ $N$_satisfies conditions (
3.1). Let’s assume that the functions _$v\in H^{2}\left(Q_{T}\right)$_and_
$v^{\ast}\in H^{2}\left(Q_{T}\right)$_are solutions of Problem 3 with the
vectors of data_
$\left(\psi_{0}\left(t\right),\psi_{1}\left(t\right)\right)$_and_
$\left(\psi_{0}^{\ast}\left(t\right),\psi_{1}^{\ast}\left(t\right)\right)$_respectively,
where_ $\psi_{0},\psi_{1},\psi_{0}^{\ast},\psi_{1}^{\ast}\in
H^{2}\left(0,T\right).$_Assume also that error estimates ( 3.5) of the
boundary data is in place. Choose an arbitrary number
_$\epsilon\in\left(0,T\right)$_. Denote_

$\lambda=\lambda\left(T,\epsilon\right)=\frac{\ln\left(T+1-\epsilon\right)}{\ln\left(T+1\right)}\in\left(0,1\right).$
(3.13)

_Then there exists a sufficiently small number_
$\nu_{0}=\nu_{0}\left(T_{1},b_{0},b_{1}\right)\in\left(0,1\right)$_and a
constant_ $C_{1}=C_{1}\left(T,b_{0},b_{1},\epsilon\right)>0,$ _both depending
only on listed parameters, such that the following stability estimate holds
for all_ $\nu\in\left(0,\nu_{0}\right):$

$\left\|v_{x}-v_{x}^{\ast}\right\|_{L_{2}\left(Q_{T-\epsilon}\right)}+\left\|v-v^{\ast}\right\|_{L_{2}\left(Q_{T-\epsilon}\right)}\leq$
(3.14) $\leq
C_{1}\left(1+\left\|v-v^{\ast}\right\|_{H^{2}\left(Q_{T}\right)}\right)\exp\left[-\left(\ln\nu^{-1/2}\right)^{\lambda}\right].$

Below $C=C\left(T,b_{0},b_{1}\right)>0$ and
$C_{1}=C_{1}\left(T,a_{0},b_{1}\right)>0$ denote different constants depending
only on listed parameters.

Corollary 4.1 (uniqueness). _Let the coefficient_ $b\left(x,t\right)$_of the
operator_ $N$_satisfies conditions ( 3.1). Then Problem 3 has at most one
solution (uniqueness)_.

Proof.__ If $\nu=0,$ then (3.14) implies that
$v\left(x,t\right)=v^{\ast}\left(x,t\right)$ in $Q_{T-\epsilon}.$ Since
$\epsilon\in\left(0,T\right)$ is an arbitrary number, then
$v\left(x,t\right)\equiv v^{\ast}\left(x,t\right)$ in $Q_{T}.$ $\square$

Theorem 4.3 (existence and uniqueness of the minimizer). _Let functions_
$\psi_{0}\left(t\right),\psi_{1}\left(t\right)\in H^{2}\left(0,T\right).$_Let_
$W$_be the set defined in ( 3.7). Then there exists unique minimizer
_$v_{\min}\in W$_of functional ( 3.6) and_

$\left\|v_{\min}\right\|_{H^{2}\left(Q_{T}\right)}\leq\frac{C}{\sqrt{\beta}}\left(\left\|\psi_{0}\right\|_{H^{2}\left(0,T\right)}+\left\|\psi_{1}\right\|_{H^{2}\left(0,T\right)}\right).$
(3.15)

In the theory of Ill-Posed Problems, this minimizer $v_{\min}$ is called
“regularized solution” of Problem 3 [9]. According to the theory of Ill-Posed
problems, it is important to establish convergence rate of regularized
solutions to the exact one $v^{\ast}.$ In doing so, one should always choose a
dependence of the regularization parameter $\beta$ on the noise level $\nu,$
i.e. $\beta=\beta\left(\nu\right)\in\left(0,1\right)$ [9].

Theorem 4.4 (convergence rate of regularized solutions). _Let_ $v^{\ast}\in
H^{2}\left(Q_{T}\right)$_be the solution of Problem 3 with the noiseless data_
$\left(\psi_{0}^{\ast}\left(t\right),\psi_{1}^{\ast}\left(t\right)\right).$_Let
functions_ $\psi_{0},\psi_{1},\psi_{0}^{\ast},\psi_{1}^{\ast}\in
H^{2}\left(0,T\right).$_Let_ $v_{\min}\in W$_be the unique minimizer of
functional ( 3.6) on the set _$W$_. Assume that error estimates ( 3.5) hold.
Choose an arbitrary number _$\epsilon\in\left(0,T\right)$_. Let_
$\lambda=\lambda\left(T,\epsilon\right)\in\left(0,1\right)$ _be the number
defined in ( 3.13) and let_

$\beta=\beta\left(\nu\right)=\nu^{2},$ (3.16)

_Then there exists a sufficiently small number_
$\nu_{0}=\nu_{0}\left(T,b_{0},b_{1}\right)\in\left(0,1\right)$_depending only
on listed parameters such that the following convergence rate of regularized
solutions_ $v_{\min}$ _holds for all_ $\nu\in\left(0,\nu_{0}\right):$__

$\left\|\partial_{x}v_{\min}-\partial_{x}v^{\ast}\right\|_{L_{2}\left(Q_{T-\epsilon}\right)}+\left\|v_{\min}-v^{\ast}\right\|_{L_{2}\left(Q_{T-\epsilon}\right)}$
(3.17) $\leq
C_{1}\left(1+\left\|v^{\ast}\right\|_{H^{2}\left(Q_{T}\right)}+\left\|\psi_{0}^{\ast}\right\|_{H^{2}\left(0,T\right)}+\left\|\psi_{1}^{\ast}\right\|_{H^{2}\left(0,T\right)}\right)\exp\left[-\left(\ln\nu^{-1/2}\right)^{\lambda}\right].$

###  3.3 Trading Strategy:

We use minimizers obtained from the method of Quasi-Reversibility to build a
strategy for trading options. Let’s define

$REAL(0)=\frac{u_{a}(0)+u_{b}(0)}{2}$ (3.18)

$REAL(\tau)=\frac{u_{a}(\tau)+u_{b}(\tau)}{2}$ (3.19)

$EST(\tau)=u_{\beta}(1/2,k\tau)$ (3.20)

or if it was not applied dimensionless

$EST(\tau)=u_{\beta}(\frac{s_{a}+s_{b}}{2},k\tau)$ (3.21)

where $k=1$

Here $EST(\tau)$ means minimizer.

Let’s buy an option if the following holds

$EST(\tau)\geq REAL(0)$ (3.22)

The predicted outcome of option trade is Positive if

$EST(\tau)\geq REAL(0)$ (3.23)

Definition 1.

It is True Positive if

$EST(\tau)\geq REAL(0)$ (3.24)

and

$REAL(\tau)\geq REAL(0)$ (3.25)

Definition 2.

It is True Negative if

$EST(\tau)<REAL(0)$ (3.26)

and

$REAL(\tau)<REAL(0)$ (3.27)

Definition 3.

It is False Positive if

$EST(\tau)\geq REAL(0)$ (3.28)

and

$REAL(\tau)<REAL(0)$ (3.29)

Definition 4.

It is False Negative if

$EST(\tau)<REAL(0)$ (3.30)

and

$REAL(\tau)\geq REAL(0)$ (3.31)

The accuracy of trading strategy is defined as

$Accuracy=\frac{TP+TN}{\sum options}$ (3.32)

where $TP$ is a summation of True Positive and $TN$ is a summation of True
Negative and $\sum options$ is a summation of options in data set.

The precision of trading strategy is defined as

$Precision=\frac{TP}{TP+FP}$ (3.33)

where $FP$ is a summation of False Positive.

The recall of trading strategy is defined as

$Recall=\frac{TP}{TP+FN}$ (3.34)

where $FN$ is a summation of False Negative.

The average relative error of trading strategy is defined as

$Error=\frac{1}{N}\sum|\frac{EST(\tau)-REAL(\tau)}{REAL(\tau)}|$ (3.35)

##  4 Application of Neural Network Machine Learning

The Black-Scholes equation gives fair value of options in perfect market.
However, real options prices contain some level of noise. We try to filter
mispredictions (i.e. where minimizers result in False Positive or False
Negative) caused by input noise using Machine Learning to improve accuracy,
precision and recall of the trading strategy. We built a neural network with
13 element input vector and 3 fully connected hidden layers. (See Fig 1).
Input vector consists of minimizers (for $t=\tau,2\tau$) obtained from the
method of Quasi-Reversibility, stock ask and bid price (for $t=0$), option ask
and bid price and volatility (for $t=-2\tau,-\tau,0$).

Figure 1:

All vectors and labels are split into three parts: training, validation and
test sets. The training set is used for weight learning. Validation set is
used for tuning of the neural network hyper-parameters. Test set is for
generating the outcomes of trading strategy.

We collected historical option and stock prices along with implied volatility
on companies consisting of Russel 2000 index [6].

Table 1.

Set | Dates | Number of options
---|---|---
Training | 2016/09/14-2018/05/31 | 132,912
Validation | 2018/06/01-2018/06/29 | 13,401
Test | 2018/07/02-2018/08/17 | 23,549

We compared the profitability of the trading strategy based on the original
minimizer set with the profitability of the output of Machine Learning.

###  4.1 Machine Learning Input Vector Normalization

$\mu=\frac{u_{a}(0)+u_{a}(-\tau)+u_{a}(-2\tau)+u_{b}(0)+u_{b}(-\tau)+u_{b}(-2\tau)}{6}$
(4.1)

$op_{n}=\frac{u(t)-\mu}{\sigma}$ (4.2)

$s_{n}=\frac{(s-st)-\mu}{\sigma}$ (4.3)

where $op_{n}$ is a normalized option price, $s_{n}$ is a normalized stock
price normalization, $s$ is the stock, $st$ is the strike and $\sigma$ is the
standard deviation.

###  4.2 Binary classification

Supervised Machine Learning has been applied to the neural network for the
Cross Entropy Loss function with regularization:

$L(\theta)=\frac{1}{m}\sum_{i=1}^{m}[-y^{(i)}\log(h_{\theta}(x^{(i)})-(1-y^{(i)})\log(1-h_{\theta}(x^{(i)}))]+\frac{\lambda}{2m}\sum_{j=1}^{n}\theta_{j}^{2}$
(4.4)

Where $\theta$ are weights which are optimized by minimizing the loss function
using the method of gradient descent. $\lambda$ is a parameter of
regularization. $x^{(i)}$ is our normalized 13 - dimensional vectors.
$h_{\theta}$ is output of the neural network. $m$ is the number of vectors in
the training set. $y^{(i)}$ is our labels (the ground truth). The trading
strategy is defined by

$H_{c}=\begin{cases}1,&\text{if }h_{\theta}>c\\\
0,&\text{otherwise}.\end{cases}$ (4.5)

where $c$ is the threshold obtained by maximizing accuracy on validation set.
The labels are set to 1 for profitable trades and 0 otherwise.

###  4.3 Regression model

Similarly, instead of using binary classification, we can use the same ML
architecture to predict the option price for tomorrow ($\tau$). We have the
same input features as the classification neural network. Regression learning
uses mean squared error as the loss function:

$L(\bar{h}_{\theta},\bar{y})=\frac{1}{m}\sum_{n}(\bar{h}_{\theta}-\bar{y}_{n})^{2}$
(4.6)

Where $m$ is the size of the data set, $\bar{h}_{\theta}$ is the predicted
value and $\bar{y}$ is the real value ($Real(\tau)$).

##  5 Results

The following graph shows the accuracy of the results on validation set. We
use it to determine the optimal value of hyper-parameter $c$ (the threshold
value of binary classification).

Figure 2: Accuracy. Threshold dependency

Observation 1.

The accuracy was improved by both Machine Learning methods compared to the
method of Quasi-Reversibility. Based on this graph we set $c=0.5.$

The next graph presents Recall and Precision diagram built on validation data
set.

Figure 3: Precision vs Recall. Threshold dependency. Here X indicates the
position of Recall and Precision produced by Regression Neural Network and
$\Diamond$ indicates the position of Recall and Precision produced by QRM.

Observation 2.

Binary Classification and Regression produced similar results that improved
both precision and recall compared to the method of Quasi-Reversibility.

Further we divided our test data into bins where (horizontal axis, see Figure
4) each bin is determined by $\frac{s-st}{s}$ with step size $0.1.$ for each
bin we calculated precision (see Figure 4).

Figure 4: Binary Classification, Regression NN and the method of Quasi-
Reversibility.

Observation 3.

To our surprise, when stock price was close to the strike price Machine
Learning and the method of Quasi-Reversibility give similar precision (bin
$0$). With stock price diverging from the strike price Machine Learning
produced better precision.

The following tables summarize the accuracy, precision and recall for all
methods on test data.

Table 2. Final results on Test Data.

Method | Accuracy | Precision | Recall | Error
---|---|---|---|---
QRM | 49.77% | 55.77% | 52.43% | 12 %
Binary Classification | 56.36% | 59.56% | 70.22% | NA
Regression NN | 55.42% | 60.32% | 61.29% | NA

Table 3. Percentages of options with profits/losses for three different
methods.

Method | Profitable options | Options with loss
---|---|---
QRM | 55.77% | 44.23 %
Binary Classification | 59.56% | 40.44%
Regression NN | 60.32% | 39.68%

##  6 Summary

To predict prices of stock options, we used two empirical mathematical models
for Black-Scholes equation. The results achieved by solving the equation
forwards in time (as an ill-posed problem) and applying Supervised Machine
Learning (Binary Classification and Regression Neural Network, and using these
methods with the real market data, show that this methodology produce
promising results, potential applications within real-world trading and
investment strategies.

The comparison of our methods resulted in the following two conclusions:

  1. 1.

The predictions of the method of Quasi-Reversibility ended up being profitable
for $55.77\%$ of the options. Compare this to a $59.56\%$ profitability rate
for the Binary Classification method, and a $60.32\%$ profitability rate for
Regression Neural Network, used on the same data set and with the same trading
strategy.

  2. 2.

As shown by figures in section 5, option price forecasting using Machine
Learning gives us significant accuracy and profit improvements over the method
of Quasi-Reversibility. However, when stock price is close to the strike price
both models give similar results.

The authors hypothesize that options traders can generate significant profits
using trading strategies reliant on predictions generated with these methods.

## References

  * [1] M.V. Klibanov A.B. Bakushinskii and N.A. Koshev, _Carleman weight functions for a globally convergent numerical method for ill-posed Cauchy problems for some quasi- linear PDEs_ , Nonlinear Analysis: Real World Applications, 34:201–224, 2017.
  * [2] https://bloomberg.com.
  * [3] M. V. Klibanov and A. G. Yagola, _Convergent numerical methods for parabolic equations with reversed time via a new Carleman estimate_ , Inverse Problems, 35, 2019.
  * [4] M.V. Klibanov, _Carleman estimates for the regularization of ill-posed Cauchy problems_ , Appl. Numer. Math., 94, 46-74, 2015\.
  * [5] M.V. Klibanov A.V. Kuzhuget and K.V. Golubnichiy, _An ill-posed problem for the Black-Scholes equation for a profitable forecast of prices of stock options on real market data_ , Inverse Problems, 32(1), 2016.
  * [6] https://money.cnn.com/data/markets/russell.
  * [7] S. E. Shreve, _Stochastic Calculus for Finance II. Continuous - Time Models_ , Springer, 2003
  * [8] M.M. Lavrent’ev, V.G. Romanov and S.P. Shishatskii, _Ill-Posed Problems of Mathematical Physics and Analysis_ , Providence, RI: American Mathematical Society, 1986.
  * [9] A. N. Tikhonov, A. V. Goncharsky, V. V. Stepanov and A. G. Yagola, _Numerical Methods for the Solution of Ill-Posed Problems_ , Kluwer Academic Publishers Group, Dordrecht, 1995.
  * [10] M.V. Klibanov, A. A. Shananin, K. V. Golubnichiy and S. M. Kravchenko, _Forecasting Stock Options Prices via the Solution of an Ill-Posed Problem for the Black-Scholes Equation_ , arXiv preprint arXiv:2202.07174.
  * [11] T. Bjork, _Arbitrage Theory in Continuous Time_ , Oxford University Press, 1999.


# Optical-coherence-tomography-based deep-learning scatterer-density estimator
using physically accurate noise model

Thitiya Seesan  1 Pradipta Mukherjee  1 Ibrahim Abd El-Sadek  1,2 Yiheng Lim
1 Lida Zhu  1 Shuichi Makita  1 and Yoshiaki Yasuno1  1Computational Optics
Group, University of Tsukuba, Tsukuba, Ibaraki 305-8573, Japan
2Department of Physics, Faculty of Science, Damietta University, New Damietta
City 34517, Damietta, Egypt<EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

https://cog-news.blogspot.com/

We demonstrate a deep-learning-based scatterer density estimator (SDE) that
processes local speckle patterns of optical coherence tomography (OCT) images
and estimates the scatterer density behind each speckle pattern. The SDE is
trained using large quantities of numerically simulated OCT images and their
associated scatterer densities. The numerical simulation uses a noise model
that incorporates the spatial properties of three types of noise, i.e., shot
noise, relative-intensity noise, and non-optical noise. The SDE’s performance
was evaluated numerically and experimentally using two types of scattering
phantom and in vitro tumor spheroids. The results confirmed that the SDE
estimates scatterer densities accurately. The estimation accuracy improved
significantly when compared with our previous deep-learning-based SDE, which
was trained using numerical speckle patterns generated from a noise model that
did not account for the spatial properties of noise.

##  1 Introduction

Optical coherence tomography (OCT) is a non-invasive imaging modality that is
used to provide high-resolution structural information about biological
tissues[1, 2, 3]. The anatomical images provided by OCT are used widely in
clinical diagnosis [4].

In addition to anatomical investigations, OCT-based assessments of the optical
properties of tissue have also been studied and have provided useful
biomarkers. One commonly investigated optical property of biological tissues
is the attenuation coefficient (AC) [5, 6], and the AC is considered to be
related to the tissue density. AC measurements are useful in a wide variety of
applications, including investigation of tumor spheroid necrosis [7, 8] and
distinguishing between normal and cancerous tissues [9, 10, 11]. However, AC
measurements are strongly influenced by the measurement configuration and
conditions, which including the system confocality, the depth position of the
focus, and the presence of aberrations [12, 13, 14], and these factors can
limit the accuracy and reliability of the AC measurements. Although several
methods have been proposed to compensate for these effects, they generally
require hard-wired assumptions to be made [12] or multiple measurements to be
performed [15].

Rather than use the AC to assess biological tissue density, we have proposed a
deep-learning based method that estimates the scatterer density of the tissue
directly [16, 17, 18]. In this work, we denote this method as the scatterer
density estimator (SDE). The SDE analyzes the local spatial patterns of an OCT
intensity image, i.e., the speckle pattern, using a convolutional neural
network (CNN) and then estimates the scatterer density. The CNN was trained
using fully numerically simulated OCT speckle patterns that were generated by
a simple scalar-optics-based OCT simulator. This approach provides significant
amounts of training data, and the data sets reflect the variety of the
parameters involved in OCT imaging, including the resolution, the signal-to-
noise ratio (SNR), and the sample scatterer density.

In our previous study[18], we noted that the CNN can estimate two primary
quantities: the speckle contrast and the resolution. Hillman et al.showed that
the local speckle contrast of OCT has a monotonic and negative relationship
with the number of scatterers contained with in a three-dimensional (3D)
resolution volume, which is designated the effective number of scatterers
(ENS) [19]. In addition, Kurokawa et al.demonstrated that the axial and
lateral resolutions of OCT can be estimated based on a local speckle pattern
[20]. Therefore, we have anticipated that the CNN will primarily estimate the
speckle contrast and the resolutions, and have then estimated the scatterer
density based on these quantities.

This previous SDE has provided promising results, but its scatterer density
estimation accuracy was limited when it was applied to real experimental data,
particularly if the scatterer density or the SNR is low [18]. By considering
the possible mechanism for the estimator described above, we hypothesized that
the low estimation accuracy for both the speckle contrasts and the resolutions
may be limiting factors for the scatterer density. In addition, the estimation
accuracy of the scatterer density falls when the SNR is low, i.e., in cases in
which the noise forms a significant component of the OCT image. As a result,
we hypothesized that our simulation-generated training data (i.e., the speckle
patterns) did not actually reflect the physical properties of the noise
accurately. This inaccuracy in the noise modeling process may have then caused
the inaccuracies in both the speckle contrast and the resolution estimations.

In this paper, we introduce a new physically realistic noise model and
demonstrate improved accuracy in the SDE. This new noise model incorporates
the spatial properties of the OCT noises including shot noise, relative
intensity noise (RIN), and non-optical noise. The accuracy of the SDE is then
validated using numerically generated OCT images, along with OCT images of two
types of scattering phantoms and in vitro tumor spheroids. Our results
demonstrate significant improvement in the scatterer density estimation
accuracy, particularly under low SNR conditions.

##  2 Principle

###  2.1 Neural network-based scatterer density estimator

Our CNN-based SDE estimates the scatterer density from a small spatial pattern
of an OCT image, i.e., from a local speckle pattern with an image size of 16
$\times$ 16 pixels. To train the CNN model, we used numerical speckle patterns
that were generated by a simple scalar-optics-based OCT simulator, which is
designated a “speckle generator.” The speckle generator generates a speckle
pattern from an arbitrary parameter set that consists of the scatterer
density, the axial and lateral resolutions, and the SNR. The speckle pattern
is generated by convolving a 3D scatterer distribution map with a 3D complex
point spread function and subsequently adding complex Gaussian noise. The
details of the speckle generator are presented in Section 3.1.

In our previous study, the noise was assumed to be fully spatially
decorrelated. In other words, the noise signals at the individual pixels are
assumed to be fully independent of each other. In the new SDE, we introduce a
more physically accurate noise model that accounts for the spatial correlation
properties of the noise, as will be described in the next section.

###  2.2 Noise Model

Similar to our previous noise model, the new noise model has a complex
Gaussian distribution, i.e., both the real and imaginary parts of the noise
are normally distributed with a zero mean. However, our new noise model
accounts for the different spatial correlation properties of the three types
of OCT noise, comprising shot noise, RIN, and non-optical noise. These three
noise types have different spatial extents, and these differences originated
from the different dependencies of these noise types on the light source
spectrum.

When the light source spectral intensity is denoted by $S(k)$, the noise
within the wavenumber domain can be described as

$N_{a}(k)=c_{1}\sqrt{S(k)}N_{\mathrm{sh}}(k)+c_{2}S(k)N_{\mathrm{RIN}}(k)+c_{3}N_{\mathrm{no}}(k),$
(1)

where $k$ is the wavenumber, and $N_{\mathrm{sh}}$, $N_{\mathrm{RIN}}$, and
$N_{\mathrm{no}}$ are normally distributed random variables that correspond to
the shot noise, the RIN, and the non-optical noise, respectively. All these
variables have a zero mean and the same standard deviations, but they are
independent of each other. In addition, the constants $c1$, $c2$, and $c3$ are
proportionality constants.

As the equation shows, the amplitude expectation of the shot noise (i.e., the
first term on the right-hand side) is a function of $k$, and is proportional
to the square root of the spectral intensity at each value of $k$. This is
because the shot noise energy, which is the squared power of the amplitude, is
proportional to the spectral intensity. Similarly, the amplitude expectation
of the RIN (i.e., the second term) is also a function of $k$, and is also
proportional to the spectral intensity. This corresponds to the fact that the
RIN energy is proportional to the squared power of the spectral intensity. In
contrast, the non-optical noise (i.e., the third term) is independent of the
light source spectrum and its amplitude expectation remains constant over the
range of $k$.

The spectral-intensity dependence of the noise causes a spatial correlation of
the noise along the depth direction. The noise in the depth domain can be
determined by taking the Fourier transform of Eq. (1) as follows

$\mathcal{F}^{-1}\left[{N_{a}(z)}\right]=c_{1}\mathcal{F}^{-1}\left[{\sqrt{S(k)}}\right](z)*N_{\mathrm{sh}}^{\prime}(z)+c_{2}\mathcal{F}^{-1}\left[{S(k)}\right](z)*N_{\mathrm{RIN}}^{\prime}(z)+c_{3}N_{\mathrm{no}}^{\prime}(z)$
(2)

where $z$ is the depth and is a Fourier pair of $k$.
$\mathcal{F}^{-1}\left[{\quad}\right]$ denotes an inverse Fourier transform
and $N_{\mathrm{sh}}^{\prime}(z)$, $N_{\mathrm{RIN}}^{\prime}(z)$, and
$N_{\mathrm{no}}^{\prime}(z)$ are the Fourier transforms of
$N_{\mathrm{sh}}(k)$, $N_{\mathrm{RIN}}(k)$, and $N_{\mathrm{no}}(k)$,
respectively. Because $N_{\mathrm{sh}}(k)$, $N_{\mathrm{RIN}}(k)$, and
$N_{\mathrm{no}}(k)$ are all normally distributed Gaussian noise,
$N_{\mathrm{sh}}^{\prime}(z)$, $N_{\mathrm{RIN}}^{\prime}(z)$, and
$N_{\mathrm{no}}^{\prime}(z)$ become normally distributed complex Gaussian
noise. As the equation shows, the shot noise is convolved with the Fourier
transform of $\sqrt{S(k)}$, and the result is the spatial correlation of the
shot noise along the depth direction. Similarly, the RIN is convolved with the
Fourier transform of $S(k)$ and this operation also results in a spatial
correlation along the depth direction.

Notably, because $\sqrt{S(k)}$ can be wider than $S(k)$,
$\mathcal{F}^{-1}\left[{\sqrt{S(k)}}\right](z)$ can then be narrower than
$\mathcal{F}^{-1}\left[{S(k)}\right](z)$. This suggests that the spatial
correlation distance of the shot noise is shorter than that for the RIN. In
contrast to the shot noise and the RIN, i.e., the optical noise types, the
non-optical noise (the third term on the right-hand side) shows no spatial
correlation. The spatial correlation distances for the three types of noise
contained in the new noise model and in the old noise model used in our
previous SDE [18] are summarized in Table 1, in addition to the correlation
distance of the OCT signal.

Table 1:  Summary of spatial correlation distance characteristics of three
noise types and the OCT signal. In the new noise model, each type of noise has
different spatial correlation properties, whereas these properties [18] are
identical for all noise types in the old noise model.

In this study, we consider scanning OCT, rather than a full-field OCT. Because
each A-line of the scanning OCT is acquired at a different time, none of the
three noise types have spatial correlation along the lateral direction. We
expect that the different spatial properties of the noises and the OCT signal
will be used by the CNN to estimate the speckle contrast and the resolutions
accurately. This will then enable more accurate estimation of the scatterer
density than our old CNN model, as will be shown in Section 4.

##  3 Implementation and validation method

###  3.1 Speckle generation and CNN-model training

The first step in the CNN model training process is the generation of the OCT
speckle patterns, i.e., the training datasets. This step is identical to that
described in Section 2.1 of [18], with the exception of the noise model.
Because the details of this numerical speckle pattern generation process have
been described elsewhere, we describe the process only briefly here. In this
process, we first generate the 3D numerical fields using randomly distributed
scatterers. Here, the field has a size of 128 $\times$ 128 $\times$ 128 pixels
(31.2 $\mu$m $\times$ 31.2 $\mu$m for the lateral directions and 115.8 $\mu$m
for the axial direction), and the pixels with a scatterer have an amplitude of
unity but a random phase. This random phase represents the sub-wavelength
depth position of the scatterer. The pixels without a scatterer have zero
values. The 3D field of scatterers is then convolved with a 3D complex point
spread function of the OCT, which is assumed to have a 3D Gaussian
distribution and its two lateral resolutions are assumed to be identical. The
3D numerical field is then down-sampled from 128 $\times$ 128 $\times$ 128
pixels to 16 $\times$ 16 $\times$ 16 pixels, in keeping with the original
physical field size. The pixel size after down-sampling is then 1.95 $\mu$m
for the lateral directions and 7.24 $\mu$m for the axial direction. After
down-sampling of the field, we added the shot noise, the RIN, and the non-
optical noise because it follows from the noise model described in Section
2.2.

To train the CNN model, we generated 80,000 3D speckle patterns and then
extracted 1,280,000 2D cross-sectional speckle patterns with dimensions of 16
$\times$ 16 pixels. Each 3D speckle pattern was generated with different and
randomly selected values for the scatterer density, the resolutions, and the
noise energies of the shot noise, the RIN, and the non-optical noise. The
axial and lateral resolutions here are independent of each other and range
from 3 to 30 $\mu$m. The noise energies for each type of noise are randomly
selected as the SNR with respect to each noise type in the range from 0 to 100
dB. The scatterer densities range from 0 to 0.2387 scatterers/$\mu$m3 (0 to
6.5723 scatterers/pixel).

The CNN consists of three convolutional and max-pooling layer pairs followed
by two fully connected layers; this structure is identical to that described
in Section 2.2.2 of Ref. [18]. The model was trained to minimize the mean
squared error (MSE) between the ground-truth scatterer density values and the
network outputs by using the Adam optimizer [21] at a learning rate of 10-4.
The batch size of the training was set at 32. Before being input into the
model, the input image intensity was normalized into a [0, 1] range. The
validation set used for the training consists of 100 2D speckle patterns that
were extracted from 100 3D speckle patterns; these patterns were generated
independently from the training dataset.

###  3.2 Evaluation method

The performances of the SDE were validated bothnumerically and experimentally.
The numerical validation used numerically generated OCT speckle patterns,
whereas the experimental validations used two types of scattering phantom
(comprising Intralipid phantoms and microsphere phantoms) and in vitro tumor
spheroid samples. For comparison, we also evaluated our previous SDE [18],
which is identical to the new SDE with the exception of the noise model used
to generate the training dataset. The process details are described in the
following sections.

####  3.2.1 Numerical validation

For the numerical validation, we used 100 2D OCT speckle patterns that were
generated using the same method that was used to generate the training
dataset. Each of these 100 2D speckle patterns was extracted from different 3D
speckle patterns. Therefore, all the 2D speckle patterns are independent of
each other and are based on different resolutions, SNR values, and scatterer
densities. Note that the noise model used for this generation process is the
physically accurate noise model that was described in Section 2.2. Because we
know the true density, we can determine the estimation accuracy directly via
this numerical validation process.

####  3.2.2 Validation using scattering phantoms

For the phantom-based experimental validations, we used two scattering phantom
types, i.e., Intralipid phantoms and microsphere phantoms.

The Intralipid phantoms are composed of Intralipid solutions with various
concentrations of 1%, 2%, 4%, 6%, 8%, and 10% (v/v concentration). These
Intralipid solutions are made from a 20% stock solution of Intralipid (IL-20,
Sigma-Aldrich I141). Three phantoms were fabricated for each concentration,
giving a total of 18 phantoms. Note that these phantoms are identical to the
phantoms used in our previous study [18], and the same measurement datasets
that were used in the previous study were also reused in this study.

The microsphere phantoms are formed using a mixture of polystyrene
microspheres and agarose gel. When compared with the Intralipid droplets, the
polystyrene microspheres have well controlled and known sizes and refractive
index values, and thus these microspheres are frequently used as phantoms in
the biomedical optics field [22]. In our study, the microsphere-based phantom
is particularly useful because we can compute the true scatterer densities
from the product-specific particle concentrations and diameter in each case.

Six phantoms with volume concentrations of 0.1%, 0.2%, 0.4%, 0.6%, 0.8%, and
1.0% were prepared by mixing 1-$\mu$m diameter polystyrene microspheres
(89904-10ML-F, Sigma-Aldrich) with the agarose gel (A1296-100G, Sigma-
Aldrich). These concentrations corresponded to scatterer densities of 0.0191,
0.0382, 0.0764, 0.1146, 0.1528, and 0.1910 scatterers/$\mu$m3, respectively.
The scatterer density was computed here using

$\mathrm{Scatterer\,\,\,density}=\frac{6\sigma}{\pi d^{3}},$ (3)

where $\sigma$ is the volume concentration of the scattering medium and $d$ is
the microsphere diameter. The mixture was poured into an acrylic container
with a thickness of 1 mm. The container was then refrigerated for 1 h at 5 ℃.
Note that these phantoms were fabricated by following the protocol described
in [23]. Three phantoms were made with each concentration, and thus 18
phantoms were fabricated in total.

Both the Intralipid phantoms and the microsphere phantoms were measured using
a swept-source OCT with a probe beam wavelength of 1.3 $\mu$m. (The OCT system
is described in detail in Section 3.2.5.) For the measurements, we
intentionally attenuated the probe beam by applying a variable neutral density
(ND) filter to alter the SNR. Each phantom was measured with three different
round-trip attenuations, which were 0, -5.4, and -11.8 dB for the Intralipid
phantoms, and 0, -5.6, and 12.0 dB for the microsphere phantoms. The probe
power on the sample was 12 mW with 0-dB attenuation.

####  3.2.3 Human breast cancer spheroid

Tumor spheroids composed of human breast adenocarcinoma (MCF-7 cell line) were
used to evaluate the SDEs from a biomedical application perspective. The cells
were cultured for 15 days and spheroids with a size of approximately 500
$\mu$m were formed. The scatterers within the cells are believed to consist of
cell nuclei and organelles. For the measurements, each spheroid was extracted
from the culturing environment and then placed in a room-temperature culture
medium without CO2 supply.

We performed two types of measurements. The first type involves longitudinal
hourly measurements for up to 28 h, where the measurements were performed
using a 2D cross-sectional scan protocol. Each cross-sectional image (i.e.,
B-scan) consists of 512 A-lines. Note that this experiment was performed
initially for the study described in Ref. [24], and that the same data set was
used in our previous SDE study [18].

The second type is a 3D measurement taken at two longitudinal time points of 0
h and 20 h. A 3D datasets consists of 512 $\times$ 128 A-lines. This
measurement was originally performed for the study detailed in Ref. [25].

####  3.2.4 Statistical analysis

For the numerical validation and the polystyrene microsphere phantom
experiment, the agreement between the estimates and the ground truth was
evaluated statistically via intraclass correlation (ICC). Here, a higher
(i.e., closer to 1.0) ICC indicates a better estimate. This statistical
analysis was performed using the intraclass_corr( ) function in the statistics
library of Python (Pingouin 0.5.3) in Python 3.7.

####  3.2.5 OCT devices and measurement protocol

A polarization-sensitive Jones matrix swept-source OCT system [26, 27] was
used to perform the experimental validations. The probe wavelength was 1.3
$\mu$m and the measurement speed was 50,000 A-lines/s. The axial resolution
and the pixel separation were 14.1 $\mu$m and 7.24 $\mu$m (both in the
tissue), respectively, while the lateral resolution and the pixel separation
were 18.1 $\mu$m and 1.95 $\mu$m, respectively. Note that the pixel
separations were identical to those of the numerically generated speckle
patterns. The OCT image used in this study is a coherent composition of
multiple polarization channels that is nearly identical to a standard non-
polarization-sensitive OCT image.

##  4 Results

###  4.1 Numerical validation

Figure 1: Numerical validation of the new SDE (red circles) and the old SDE
(blue $\times$). The horizontal and vertical axes correspond to the set
scatterer density (the ground truth) and the estimated densities,
respectively. The black solid lines represents the perfect estimate. The new
SDE gives reasonable estimates, whereas the estimates from the old SDE are
downshifted significantly from the ground truth.

Figure 1 shows the numerical validation results, where the estimated scatterer
densities are plotted versus the set (ground truth) scatterer densities. The
SDE that was trained using the physically accurate noise model (the new SDE,
red circles) shows high consistency between the set and estimated scatterer
densities. In contrast, the previous SDE [18], which was trained using a
spatially-uncorrelated noise model (the old SDE, blue triangles), shows higher
variation among the estimates, and the estimated values are also downshifted
significantly from the ground truth. The ICCs between the estimates and the
ground truth were computed to be 0.975 for the new SDE and 0.722 for the old
SDE. The higher ICC of the new SDE is a quantitative demonstration of the
superior performance of the new SDE.

###  4.2 Intralipid phantom

Figure 2:  Scatterer density estimation results for Intralipid phantoms with
several volume concentrations and several probe-beam attenuations. The circles
and the triangles represent the estimates obtained when using the new and old
SDEs, respectively. The colors represent the probe beam attenuations that
correspond to the SNR of the OCT image. The error bar shows the standard
deviation for three measurements of the three phantoms. Both SDEs show high
repeatability (i.e., a small standard deviation in each case) and consistent
estimates among the different probe-beam attenuations. However, the old SDE
shows a significantly large intercept that indicates the low fidelity of the
old SDE, particularly at low concentrations. The lines are linear regression
lines computed from the data.

The results of the Intralipid phantom validation are summarized in Fig. 2. In
this figure, the estimated scatterer densities are plotted versus the volume
concentrations. The circles and triangles represent the estimates from the new
and old SDEs, respectively. The colors of the plots indicate the probe-beam
attenuations in each case. Specifically, each color corresponds to a different
SNR. The error bars represent the standard deviations among the three
measurements of the three phantoms and the lines are linear regression lines
that were computed from the data.

Both SDEs showed reasonably small variations (standard deviations) among the
phantoms, and shoewd highly linear relationships between the estimates (on the
vertical axis) and the volume concentrations (on the horizontal axis). In
addition, the estimated values were not sensitive to probe attenuation.
Specifically, both SDEs were unaffected by the SNR of the OCT image. However,
the intercept of the old SDE was as high as 0.05 $\mu$m-3 at the 0% volume
concentration, whereas that of the new SDE was close to zero (0.01 $\mu$m-3).
Because the scatterer density at the 0% volume concentration might be 0
$\mu$m-3, we can conclude that the new SDE provides significantly better
estimation accuracy when compared with the old SDE.

###  4.3 Polystyrene microsphere phantom

Figure 3:  Scatterer density estimation results for the microsphere phantoms.
The vertical axis represents the estimated scatterer density and the
horizontal axis represents the theoretically predicted scatterer density as
computed from the volume concentration of the phantom using Eq. (3). The error
bars represents the standard deviations among the three estimates of the three
samples. The black line represents the ideal estimate, while the other lines
are linear-regression lines that were computed from the data. The error bars
are almost unrecognizable because of the high repeatability of the estimates.
The new SDE (circles) provides estimates that are very close to the
theoretical predictions (black line). In contrast, the estimate of the old SDE
show significant departures from the theoretical predictions. In addition, the
estimation results are unaffected by the probe-beam attenuation. Specifically,
although the attenuations are represented by different colors, the plots and
the regression lines with the different attenuations overlap significantly and
are not really distinguishable.

The estimation results for the microsphere phantoms are plotted versus the
theoretically computed scatterer density in Fig. 3. The theoretical scatterer
density was calculated from the volume concentration of the phantom using Eq.
(3). Similar to Fig. 2, the circles and the triangles represent the estimates
from the new and old SDEs, respectively, and the colors indicate the probe-
beam attenuation levels. The error bars indicate the standard deviations among
the three measurements of the three phantoms, but the error bars are nearly
invisible here because of the very high measurement repeatability (leading to
very small standard deviations). In addition, both SDEs are nearly perfectly
independent of the SNRs, as indicated by the fact that the data points for the
different SNRs are nearly perfectly overlapping each other.

The estimates from both SDEs show highly linear relationships with the
theoretically predicted scatterer density. However, the estimates from the old
SDE (triangles) show significant departures from the ideal estimation line
(i.e., the solid black line). In contrast, the new SDE (circles) shows very
close agreement with the theoretical predictions. The ICCs between the
estimateds values and the theoretical scatterer density values were very high
at 0.982 (0-dB attenuation), 0.982 (-5.6-dB attenuation), and 0.981 (-12-dB
attenuation) for the new SDE, whereas those for the old SDE were 0.467 (0-dB
attenuation), 0.467 (-5.6-dB attenuation), and 0.460 (-11.8-dB attenuation).
The higher ICCs indicate the superior performance of the new SDE.

###  4.4 Tumor spheroids

Figure 4:  Time-lapse images of a tumor spheroid. The individual rows show (a)
the OCT intensity images, and the scatterer density images obtained by the (b)
new and (c) old SDEs. Both SDEs show a reduction in the scatterer density over
time.  Figure 5:  Mean scatterer density within the spheroid region. Both the
new (red circles) and old (blue rectangles) SDEs showed a reduction in the
mean scatterer density over time. For the time points after 10 h, the old SDE
gives constant scatterer density estimates, whereas the new SDE shows a
continuous reduction in the scatterer density. As discussed in Section 5.1,
the results from the new SDE are more plausible in this case.

Figure 4 shows hourly time-lapse images of the MCF-7 spheroid. The individual
rows show the OCT B-scans, the scatterer density images obtained using the new
SDE, and those obtained using the old SDE (in order from top to bottom). The
color scale (hue) in the scatterer density images represents the estimated
scatterer density, while the brightness corresponds to the OCT intensity. Note
that the images from the old SDE are identical to those presented in Ref.
[18]. In addition, the mean scatterer density was computed within manually
segmented spheroid regions and plotted versus time, as shown in Fig. 5. In the
figure, the red circles and the blue rectangles correspond to the new and old
SDEs, respectively.

Both SDEs showed a reduction in the scatterer density over time, but
distinctive differences between the two SDEs can be seen at the later time
points after 10 h. At these time points, the new SDE showed a continuous
reduction in the scatterer density, whereas the estimates from the old SDE
remained constant. The constant estimation from the old SDE may be caused by
the relatively large error of low SDE at lower scatterer densities, as shown
in the numerical and phantom-based validations presented in Sections 4.1, 4.2,
and 4.3. Here, we remind that the low scatterer density may cause a lower OCT
signal intensity and thus lead to the higher noise dominance. Therefore, the
scatterer density estimation can be more helpful in maintaining the noise
model accuracy for the lower scatter density cases. The continuous reduction
in the scatterer density found by the new SDE is more consistent with to our
dynamic OCT findings than the results from the old SDE; this will be addressed
later in the discussion section (Section 5.1).

Figure 6:  3D time-course images of MCF-7 spheroid. (a-d) show OCT intensity
images, and (e-h) and (i-l) are scatterer density images acquired when using
the new and old estimators, respectively. The first and second rows represent
the cross-sectional and en face images, respectively. At the zero-hour time
point, the OCT intensity did not show a clear structure within the spheroid,
whereas the en face scatterer density images showed a clear low-density core
at their center. In addition, the en face scatterer density images showed
clear difference between two time points, while this stark contrast was not
observed in the en face OCT intensity images. The scale bar represents 200
$\mu$m

Figure 6 shows cross-sectional (first row) and en face (second row) images
obtained from a volumetric dataset of a tumor spheroid. The en face images are
the slices acquired around the equator of the spheroid. The left images (a-d)
are the OCT intensity images, the middle images (e-h) are the scatterer
density images acquired using the new SDE, and the right (i-l) are those
acquired when using the old SDE.

The cross-sectional images showed similar appearances to the corresponding
images in Fig. 4. However, the en face images illustrate that the scatterer
density images are more informative than the OCT intensity images. For
example, the OCT intensity at zero-hour time point [Fig. 6(c)] only shows a
homogeneous appearance. In contrast, the scatterer density image from the new
SDE [Fig. 6(g)] revealed that the center region has a lower scatterer density
than the peripheral regions. This may be an indication of the well-known
necrotic core of MCF-7 spheroids [28, 29, 25]. Notably, the core is not
visible when using the old SDE [Fig. 6(k)]. This may occur because the
estimator accuracy of the old SDE is too low for the low scatterer density in
this case.

The en face OCT intensity images at the two time points [Fig. 6(c) and (d)] do
not show clear differences, whereas the scatterer density images [Fig. 6(g)
versus (h) and Fig. 6(k) versus (l)] show clear reductions in the scatterer
density at the late time point (20 h). The reduction in the scatterer density
at the later time point is reasonable because the dominant scatterers in these
cells are the nuclei and cell organelles, and they are resolved by cell death.
Some discussions related to this point can be found in Section 5.1. It should
also be noted that the scatterer density from the old SDE at the later time
point [Fig. 6(l)] is higher than that obtained from the new SDE [Fig. 6(h)].
This difference also can be attributed to the low accuracy of the old SDE for
low scatterer densities.

Additionally, we can see that the scatterer densities of the upper and lower
parts of the spheroid are not symmetrical. Specifically, the upper hemisphere
shows a higher scatterer density than the lower hemisphere. This issue is
discussed in greater detail in Section 5.2.

##  5 Discussion

###  5.1 Time course reduction and scatterer density

Both the 2D and 3D spheroid imaging results presentrd in Section 4.4 showed
time-course reductions in the scatterer densities. As shown in Fig. 5, the new
SDE indicates a continuous reduction in the scatterer density over 28 h, while
the scatterer density estimated using the old SDE becomes constant after 10
hours.

These reductions in the scatterer density can be understood more easily when
using dynamic OCT imaging. In our previous paper, the same datasets used in
Figs. 4 and 5 were analyzed using two dynamic OCT methods, i.e., the
logarithmic intensity variance (LIV) method and the late OCT correlation decay
speed (OCDSl) method. The mean LIV and mean OCDS$\l$ values within the entire
spheroid region were found to be decreasing continuously over 28 hours. These
reductions in the dynamic OCT signals indicate the progression of necrosis
within the spheroid [25].

The necrotic cell death causes destruction of the nuclei and the cell
organelles. Because the nuclei and the organelles are the dominant scatterers
within the cell, the progression of the necrosis causes a reduction in the
scatterer density.

According to the dynamic OCT analysis, the necrosis has progressed
continuously over 28 h; this may also suggest that the reduction in the
scatterer density also may progress over this time period. Therefore, the
results from the new SDE, which showed a continuous reduction in the scatterer
density, are more plausible than the results from the old SDE, in which the
scatterer density becomes constant after 10 h.

###  5.2 Asymmetric appearance of the scatterer density in spheroids

Cross-sectional scatterer density images of the spheroid, where Fig. 4 (b) and
(c) show asymmetric appearances between the upper and lower parts of the
spheroid. This asymmetry may partially be attributed to the occurrence of
multiple scattering in the deeper region of the tissue. Because our speckle
generator does not incorporate the multiple scattering effect, the estimation
accuracy of the SDE may decrease in the deeper region.

Several possible solutions have been proposed to address this issue. One
approach is to apply one of the available multiple-scattering rejection
methods [30, 31, 32, 33]. Another possible solution is to improve the speckle
generator to allow it to account for the multiple scattering effect. Several
models and simulation methods have been proposed that can account for multiple
scattering [34, 35, 36, 37, 38]. Although these current methods are
computationally intensive and rather too slow to generate the massive amount
of training data required, future investigations of the theoretical model and
associated simulation methods may enable development of a speckle generator
that can account for multiple scattering with reasonable computation speeds.

###  5.3 Computation time of the speckle generator and the SDE

The overall training process of the CNN-based SDE can be split into two sub-
processes, i.e., speckle generation and CNN training. The generation of 80,000
3D speckle patterns with dimensions of 16 $\times$ 16 $\times$ 16 pixels takes
approximately 3 hours when using a desktop PC equipped with an Intel Core
i7-6900K CPU and a graphics processing unit (GPU; GeForce RTX 3090 Ti,
NVIDIA). The main memory of the PC is 128 GB in size, and the GPU has 10752
cores, a 1.86 GHz boost clock, and 24 GB of memory. The speckle generator is
written in Python (version 3.7) with a GPU-compatible numerical computation
library (CuPy 12.0.0). Notably, the speckle generator is more than 11 times
faster than our previous version, which did not use a GPU. Therefore, usage of
a GPU is essential when generating large training datasets.

Subsequently, training of the CNN model takes approximately 1 h, where the CNN
is written in Python 3.7 with Keras 2.3.1 based on the TensorFlow backend.
Therefore, the entire training process takes approximately 4 hours.

Scatterer density estimation of a 512 $\times$ 402 pixel cross-sectional image
using a 16 $\times$ 16 pixels sliding window took around 1 minute, and thus
the estimation process for a volumetric image that consist of 128 B-scans
takes approximately 2 hours.

###  5.4 Compatibility of the trained model and the scanning protocol

Because the speckle generator generates a speckle pattern with a specific
pixel size, each training dataset (i.e., each set of speckle patterns) is
specific to a single scan protocol. Therefore, speckle patterns should be
generated for each different scan protocol and a specific SDE model should
then be trained using these specific speckle patterns.

The other compatibility issue that must be considered is the minimum
separation required between adjacent A-lines (i.e., the lateral pixel size).
We suspect that our new SDE uses the different lateral spatial properties of
the noise and OCT signals to distinguish them, as summarized in Table 1.
Specifically, we have assumed that the noise has fully decorrelated between
the adjacent A-lines, as described in Section 2.2. However, the OCT signals
are mutually correlated along the lateral direction with the correlation
distance around the lateral resolution. The different lateral spatial
properties of the noise and the OCT may have been used by the SDE model to
discriminate the OCT signal from the noise. However, the different lateral
extents only become pronounced if the A-line separation is smaller than the
lateral resolution. Therefore, the OCT scan protocol to be used with our SDE
should have an A-line separation that is smaller than the lateral resolution.

It is also notable that, according to Eq. (2), the RIN may have the same axial
extent as the OCT signal. Thus, if the CNN model cannot use the different
lateral properties of the signal and the noise, as in the case where the
A-line separation is far greater than the lateral resolution, the CNN model
cannot discriminate the RIN from the OCT signal. This further emphasizes the
need to consider the minimum separation required for the A-lines.

##  6 Conclusion

In this study, we demonstrated a CNN-based SDE that is trained using speckle
patterns generated by a numerical speckle generator. This speckle generator
uses a physically realistic noise model that takes the spatial properties of
shot noise, RIN, and non-optical noise into account. The SDE was examined
using numerically generated OCT speckle patterns along with the experimentally
obtained OCT images of two types of scattering phantom and in vitro tumor
spheroid samples. In these examinations, the SDE’s estimation performance was
compared with that of our previous CNN-based SDE, which was trained with
speckle patterns that did not account for the spatial properties of the noise.
All the validation experiments demonstrated the superior estimation accuracy
provided by the new SDE. In particular, the validation using microsphere
phantoms demonstrated the excellent agreement of the estimated scatterer
density with the ground truth.

In future work, this new SDE can be used in a variety of biomedical
application to assess the sub-resolution structural properties of tissues and
cells.

## Disclosures

T. Seesan, P. Mukherjee, I.A. El-Sadek, Y. Lim, L. Zhu, S. Makita, Y. Yasuno:
Sky Technology (F), Nikon (F), Kao Corp. (F), Topcon (F).

## Funding

Core Research for Evolutional Science and Technology (JPMJCR2105); Japan
Society for the Promotion of Science (21H01836, 22K04962, 22KF0058); Japan
Science and Technology Agency (JPMJFS2106); China Scholarship Council
(201908130130).

## Data, Materials, and Code Availability

Data underlying the results presented in this paper are not publicly available
at this time but may be obtained from the authors upon reasonable request.

## 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, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254, 1178–1181 (1991).
  * [2] J. G. Fujimoto, C. Pitris, S. A. Boppart, and M. E. Brezinski, “Optical coherence tomography: an emerging technology for biomedical imaging and optical biopsy,” Neoplasia 2, 9–25 (2000).
  * [3] J. G. Fujimoto, “Optical coherence tomography for ultrahigh resolution in vivo imaging,” Nat. Biotechnol. 21, 1361–1367 (2003).
  * [4] J. Chen and L. Lee, “Clinical applications and new developments of optical coherence tomography: an evidence-based review,” Clin. Exp. Optom. 90, 317–335 (2007). PMID: 17697178.
  * [5] S. Liu, “Tissue characterization with depth-resolved attenuation coefficient and backscatter term in intravascular optical coherence tomography images,” J. Biomed. Opt. 22, 1 (2017).
  * [6] S. Chang and A. K. Bowden, “Review of methods and applications of attenuation coefficient measurements with optical coherence tomography,” J. Biomed. Opt. 24, 1 (2019).
  * [7] Y. Huang, S. Wang, Q. Guo, S. Kessel, I. Rubinoff, L. L.-Y. Chan, P. Li, Y. Liu, J. Qiu, and C. Zhou, “Optical coherence tomography detects necrotic regions and volumetrically quantifies multicellular tumor spheroids,” Cancer Res. 77, 6011–6020 (2017).
  * [8] Y. Huang, J. Zou, M. Badar, J. Liu, W. Shi, S. Wang, Q. Guo, X. Wang, S. Kessel, L. L.-Y. Chan, P. Li, Y. Liu, J. Qiu, and C. Zhou, “Longitudinal morphological and physiological monitoring of three-dimensional tumor spheroids using optical coherence tomography,” J. Vis. Exp. (2020).
  * [9] E. C. C. Cauberg, D. M. de Bruin, D. J. Faber, T. M. de Reijke, M. Visser, J. J. M. C. H. de la Rosette, and T. G. van Leeuwen, “Quantitative measurement of attenuation coefficients of bladder biopsies using optical coherence tomography for grading urothelial carcinoma of the bladder,” J. Biomed. Opt. 15, 066013 (2010).
  * [10] Q. Zhao, C. Zhou, H. Wei, Y. He, X. Chai, and Q. Ren, “Ex vivo determination of glucose permeability and optical attenuation coefficient in normal and adenomatous human colon tissues using spectral domain optical coherence tomography,” J. Biomed. Opt. 17, 1050041 (2012).
  * [11] B. G. Muller, D. M. de Bruin, M. J. Brandt, W. van den Bos, S. van Huystee, D. J. Faber, D. Savci, P. J. Zondervan, T. M. de Reijke, M. P. Laguna-Pes, T. G. van Leeuwen, and J. J. M. C. H. de la Rosette, “Prostate cancer diagnosis by optical coherence tomography: First results from a needle based optical platform for tissue sampling,” J. Biophotonics 9, 490–498 (2016).
  * [12] T. G. van Leeuwen, D. J. Faber, and M. C. Aalders, “Measurement of the axial point spread function in scattering media using single-mode fiber-based optical coherence tomography,” IEEE J. Sel. Top. Quantum Electron. 9, 227–233 (2003).
  * [13] S. Stefan, K.-S. Jeong, C. Polucha, N. Tapinos, S. A. Toms, and J. Lee, “Determination of confocal profile and curved focal plane for OCT mapping of the attenuation coefficient,” Biomed. Opt. Express 9, 5084 (2018).
  * [14] K. Li, W. Liang, Z. Yang, Y. Liang, and S. Wan, “Robust, accurate depth-resolved attenuation characterization in optical coherence tomography,” Biomed. Opt. Express 11, 672 (2020).
  * [15] B. Baumann, C. W. Merkle, R. A. Leitgeb, M. Augustin, A. Wartak, M. Pircher, and C. K. Hitzenberger, “Signal averaging improves signal-to-noise in OCT images: But which approach works best, and when?” Biomed. Opt. Express 10, 5755–5775 (2019).
  * [16] T. Seesan, D. Oida, K. Oikawa, P. Buranasiri, and Y. Yasuno, “Intensity-invariant scatterer density estimation for optical coherence tomography using deep convolutional neural network,” Proc. SPIE, Biomedical Imaging and Sensing Conference 2020 11521, 115210Q (2020).
  * [17] T. Seesan, D. Oida, K. Oikawa, P. Buranasiri, and Y. Yasuno, “Sample and system parameter estimation from local speckle pattern by fully numerically trained deep convolution neural network,” Proc. SPIE, Optical Coherence Tomography and Coherence Domain Optical Methods in Biomedicine XXV 11630, 1163021 (2021).
  * [18] T. Seesan, I. Abd El-Sadek, P. Mukherjee, L. Zhu, K. Oikawa, A. Miyazawa, L. T.-W. Shen, S. Matsusaka, P. Buranasiri, S. Makita, and Y. Yasuno, “Deep convolutional neural network-based scatterer density and resolution estimators in optical coherence tomography,” Biomed. Opt. Express 13, 168 (2022).
  * [19] T. R. Hillman, S. G. Adie, V. Seemann, J. J. Armstrong, S. L. Jacques, and D. D. Sampson, “Correlation of static speckle with sample properties in optical coherence tomography,” Opt. Lett. 31, 190 (2006).
  * [20] K. Kurokawa, S. Makita, Y.-J. Hong, and Y. Yasuno, “In-plane and out-of-plane tissue micro-displacement measurement by correlation coefficients of optical coherence tomography,” Opt. Lett. 40, 2153 (2015).
  * [21] D. P. Kingma and J. Ba, “Adam: a method for stochastic optimization,” arXiv:1412.6980 [cs] (2017). ArXiv: 1412.6980.
  * [22] B. W. Pogue and M. S. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt. 11, 041102 (2006).
  * [23] A. Mustari, I. Nishidate, M. A. Wares, T. Maeda, S. Kawauchi, S. Sato, M. Sato, and Y. Aizu, “Agarose-based Tissue Mimicking Optical Phantoms for Diffuse Reflectance Spectroscopy,” J. Visualized Exp. p. 57578 (2018).
  * [24] I. Abd El-Sadek, A. Miyazawa, L. Tzu-Wei Shen, S. Makita, S. Fukuda, T. Yamashita, Y. Oka, P. Mukherjee, S. Matsusaka, T. Oshika, H. Kano, and Y. Yasuno, “Optical coherence tomography-based tissue dynamics imaging for longitudinal and drug response evaluation of tumor spheroids,” Biomed. Opt. Express 11, 6231 (2020).
  * [25] I. A. El-Sadek, A. Miyazawa, L. T.-W. Shen, S. Makita, P. Mukherjee, A. Lichtenegger, S. Matsusaka, and Y. Yasuno, “Three-dimensional dynamics optical coherence tomography for tumor spheroid evaluation,” Biomed. Opt. Express 12, 6844 (2021).
  * [26] E. Li, S. Makita, Y.-J. Hong, D. Kasaragod, and Y. Yasuno, “Three-dimensional multi-contrast imaging of in vivo human skin by Jones matrix optical coherence tomography,” Biomed. Opt. Express 8, 1290 (2017).
  * [27] A. Miyazawa, S. Makita, E. Li, K. Yamazaki, M. Kobayashi, S. Sakai, and Y. Yasuno, “Polarization-sensitive optical coherence elastography,” Biomed. Opt. Express 10, 5162 (2019).
  * [28] F. Hirschhaeuser, H. Menne, C. Dittfeld, J. West, W. Mueller-Klieser, and L. A. Kunz-Schughart, “Multicellular tumor spheroids: An underestimated tool is catching up again,” J. Biotechnol. 148, 3–15 (2010). Organotypic Tissue Culture for Substance Testing.
  * [29] E. C. Costa, A. F. Moreira, D. de Melo-Diogo, V. M. Gaspar, M. P. Carvalho, and I. J. Correia, “3D tumor spheroids: an overview on the tools and techniques used for their analysis,” Biotechnol. Adv. 34, 1427–1441 (2016).
  * [30] D. Borycki, M. Hamkało, M. Nowakowski, M. Szkulmowski, and M. Wojtkowski, “Spatiotemporal optical coherence (STOC) manipulation suppresses coherent cross-talk in full-field swept-source optical coherence tomography,” Biomed. Opt. Express 10, 2032 (2019).
  * [31] E. Auksorius, D. Borycki, and M. Wojtkowski, “Crosstalk-free volumetric in vivo imaging of a human retina with Fourier-domain full-field optical coherence tomography,” Biomed. Opt. Express 10, 6390 (2019).
  * [32] E. Auksorius, D. Borycki, P. Wegrzyn, I. Žičkienė, K. Adomavičius, B. L. Sikorski, and M. Wojtkowski, “Multimode fiber as a tool to reduce cross talk in Fourier-domain full-field optical coherence tomography,” Opt. Lett. 47, 838 (2022).
  * [33] L. Zhu, S. Makita, J. Tamaoki, A. Lichtenegger, Y. Lim, Y. Zhu, M. Kobayashi, and Y. Yasuno, “Multi-focus averaging for multiple scattering suppression in optical coherence tomography,” Biomed. Opt. Express 14, 4828–4844 (2023).
  * [34] J. M. Schmitt and A. Knüttel, “Model of optical coherence tomography of heterogeneous tissue,” J. Opt. Soc. Am. A 14, 1231–1242 (1997).
  * [35] 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).
  * [36] B. Karamata, P. Lambelet, M. Laubscher, M. Leutenegger, S. Bourquin, and T. Lasser, “Multiple scattering in optical coherence tomography I Investigation and modeling,” J. Opt. Soc. Am. A 22, 1369 (2005).
  * [37] 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).
  * [38] P. R. T. Munro, “Three-dimensional full wave model of image formation in optical coherence tomography,” Opt. Express 24, 27016–27031 (2016).