# Rational approximation on cubic hypersurfaces

David McKinnon

###### Abstract.

We compute the constant of approximation for an arbitrary rational point on an
arbitrary smooth cubic hypersurface $X$ over a number field $k$, provided that
there is a $k$-rational line somewhere on $X$. In the process, we verify the
Coba conjecture for $X$.

##  1\. Introduction

Let $X$ be a smooth projective variety defined over a number field $k$. In
[McK], the author makes a conjecture that the best rational approximations to
a rational point $P$ on a smooth projective variety can be found along some
rational curve through $P$. Such a curve is called a curve of best
approximation to $P$.

This is a little vague, obviously, but – as happens so frequently in
mathematics – we must wade through some tedious definitions to make it
precise. To wit: we must say what we mean by a good approximation.

Generally speaking, an approximation to $P$ is good (in this sense) if:

  1. (a)

the approximating point is close to $P$, and

  2. (b)

the approximating point has small height.

Let’s quantify that.

To define distance, we choose a place $v$ of $k$, and use the distance on $X$
induced by $v$. That is, we embed $X$ in projective space and restrict the
standard $v$-adic distance there to $X$. (It turns out it doesn’t matter which
embedding we pick. More on that later.)

To define the height, we choose a divisor class $L$ on $X$ – probably an ample
one – and choose some (multiplicative, relative to $k$) Weil height $H_{L}$
associated to $L$.

The tricky part is measuring the balance of these two factors. To do this, we
work in aggregate over a sequence of rational points.

Let $\\{x_{i}\\}$ be a sequence of rational points that converges $v$-adically
to $P$. The $L$-approximation constant of $\\{x_{i}\\}$ is
$\alpha(\\{x_{i}\\},P,L)$, a real number whose definition is tedious enough to
put it off to the next section. Nevertheless, $\alpha(\\{x_{i}\\},P,L)$ has
the property that the larger it is, the worse the sequence $\\{x_{i}\\}$
$v$-adically $L$-approximates $P$. It is, if you like, the cost of
approximating $P$ with the $\\{x_{i}\\}$.

We define the $L$-approximation constant $\alpha(P,L)$ to be the infimum of
$\alpha(\\{x_{i}\\},P,L)$ over all nonconstant sequences $\\{x_{i}\\}$.

Thus, the lower $\alpha(\\{x_{i}\\},P,L)$ is, the better the sequence
approximates $P$. The lower $\alpha(P,L)$ is, the better $P$ can be
approximated.

Time to deliver on my promise.

######  Conjecture 1.1 (Coba conjecture).

Let $X$ be a smooth variety defined over a number field $k$, and $P\in X$ a
$k$-rational point. Let $A$ be an ample line bundle on $X$. If
$\alpha(P,A)<\infty$, then there is a curve of best $A$-approximation to $P$.
In other words, there is a curve $C$ such that

$\alpha(P,A|_{C})=\alpha(P,A)$

where the first $\alpha$ is computed on $C$, and the second is computed on
$X$.

To paraphrase the conjecture still further: The best $A$-approximations to $P$
approximate $P$ just as well as the best $A$-approximations to $P$ that happen
to lie on $C$.

This conjecture has been verified in a large number of cases, and Vojta’s
conjecture is known to imply it in a much larger set of cases. See [McK] for a
discussion of this.

In forthcoming work of Lehmann, McKinnon, and Satriano, a framework is
proposed to deduce the  (Coba conjecture). from Vojta’s Conjecture for all
smooth projective varieties. This framework uses the machinery of the Minimal
Model Program to reduce from a general smooth variety to a Mori fibre space
(though not necessarily a smooth one). Thus, it’s natural to look to Mori
fibre spaces to verify the  (Coba conjecture).. In this paper, as warned in
the title, we consider cubic hypersurfaces. In dimension three or greater,
these are all Mori fibre spaces when they’re smooth.

In particular, we verify the  (Coba conjecture). for smooth cubic
hypersurfaces that contain a $K$-rational line, in Theorem 3.3. It is known by
work of Dietmann and Wooley (Theorem 2 of [DW]) that any smooth cubic
hypersurface of dimension at least 35 defined over a number field $k$ contains
a $k$-rational line. Work of Brandes and Dietmann (Theorem 1.1 of [BD])
improves this to dimension at least 29 when $k={\mathbb{Q}}$, so it is likely
that the hypothesis in Theorem 3.3 is satisfied in a large proportion of
cases.

The structure of the paper is quite simple. Section 2 contains some tedious
prerequisites like the definition of $\alpha$. Section 3 first considers the
cases of curves and surfaces separately in Theorem 3.1, and then proves the
higher-dimensional case in Theorem 3.3.

The author thanks Julia Brandes and Mike Roth for help and advice during the
preparation of this paper.

##  2\. Preliminaries

The purpose of this section is to make precise our terminology. Sigh. Let’s
get to it.

We start with a projective variety $X$ defined over a number field $k$. For
any line bundle $L$ on $X$, we associate a multiplicative Weil height $H_{L}$.
Details about how to do this can be found in a number of standard references,
such as section 1.1 of [Se]. There is a small ambiguity in the choice of
height, but for any two heights $H_{L}$ and $H^{\prime}_{L}$ constructed in
this way, the function $H_{L}/H^{\prime}_{L}$ is a bounded, positive function
compactly supported away from zero. And this ambiguity will not affect the
definition of $\alpha$, as we will shortly see below.

Next, we fix a place $v$ of $k$. We need to define a $v$-adic distance
function on $X$. As hinted in the introduction, we do this by choosing an
embedding of $X$ in ${\mathbb{P}}^{n}$, and restrict the usual $v$-adic
distance on ${\mathbb{P}}^{n}$ to $X$. The details of this are surprisingly
unpleasant and are covered in [MR], section 2. The important points are (a)
defining the $v$-adic distance can be done, and (b) any two distances
constructed in this way are equivalent, in that
$\operatorname{dist}/\operatorname{dist}^{\prime}$ is a bounded, positive
function compactly supported away from zero. As with the heights, this will
not affect the definition of $\alpha$.

And now for $\alpha$. Unfortunately, the path to the definition of $\alpha$ is
not a pleasingly short one. We begin by defining the approximation constant
for a single sequence.

######  Definition 2.1.

Let $\\{x_{i}\\}$ be a sequence of $k$-rational points that converges
$v$-adically to $P$. We define

$A(\\{x_{i}\\},L)=\left\\{{\gamma\in{\mathbb{R}}\mid\operatorname{dist}(P,x_{i})^{\gamma}H_{L}(x_{i})\,\,\mbox{is
bounded from above}}\right\\}.$

The constant of $L$-approximation of $\\{x_{i}\\}$ to $P$ is

$\alpha(\\{x_{i}\\},P,L)=\inf A(\\{x_{i}\\},L)$

In particular, if $A(\\{x_{i}\\},L)=\emptyset$ then
$\alpha(\\{x_{i}\\},P,L)=\infty$. We call $\alpha(\\{x_{i}\\},P,L)$ the
$L$-approximation constant of $\\{x_{i}\\}$.

The idea behind $\alpha$ is that it’s the exponent you need on
$\operatorname{dist}(x_{i},P)$ to make $\operatorname{dist}(x_{i},P)^{\alpha}$
approximately equal to $H_{L}(x_{i})$. The preceding painful definition is
merely making that precise.

Also, as advertised earlier, if we modify either the height or the distance by
a bounded, positive function compactly supported away from zero, then the
value of $\alpha$ is unchanged.

That’s $\alpha$ for a sequence. Now we need $\alpha$ for a point, independent
of the sequence. It is, unsurprisingly, the infimum of all the approximation
constants of sequences.

######  Definition 2.2.

The $L$-approximation constant of $P$ is defined to be

$\alpha(P,L)=\inf_{\\{x_{i}\\}}\alpha(\\{x_{i}\\},P,L)$

See? Told you.

Next up, a couple of results from [MR2] that will come in handy later. First,
a result about how to compute $\alpha$ on a curve.

######  Theorem 2.3 (Theorem 2.6 from [MR2]).

Let $C$ be any curve birational over $k$ to ${\mathbb{P}}^{1}_{k}$, and
$\varphi\colon{\mathbb{P}}^{1}\rightarrow C$ the normalization map. Then for
any ample line bundle $L$ on $C$, and any $P\in C(\overline{\mathbb{Q}})$ we
have the equality:

$\alpha(P,L|_{C})=\min_{q\in\varphi^{-1}(P)}d/r_{q}m_{q}$

where $d=\deg(L)$, $m_{q}$ is the multiplicity of the branch of $C$ through
$P$ corresponding to $q$, and

$r_{q}=\begin{cases}0&\text{if $\kappa(q)\not\subseteq k_{v}$}\\\ 1&\text{if
$\kappa(q)=k$}\\\ 2&\text{otherwise.}\end{cases}$

where $\kappa(q)$ denotes the field of definition of $q$.

The other result is a Liouville-type result that provides lower bounds for
$\alpha$. This is a special case of Theorem 3.3 from [MR2].

######  Theorem 2.4.

Let $X$ be an algebraic variety defined over $k$, $P\in X(k)$ any $k$-rational
point.

Let $\pi\colon\tilde{X}\to X$ be the blowup of $X$ at $P$ with exceptional
divisor $E$. Let $L$ be a nef line bundle on $X$, and $\gamma>0$ a rational
number such that $L_{\gamma}:=\pi^{*}L-\gamma E$ is in the effective cone of
$\tilde{X}$. Finally let $B^{\prime}$ be the stable base locus of $L_{\gamma}$
and set $B=\pi(B^{\prime})$.

  1. (a)

There is a positive real constant $M$ such that for all $y\in X(k)-B(k)$, we
have $H_{L}(y)\mbox{dist}_{v}(x,y)^{\gamma}\geq M$.

  2. (b)

For any sequence $\\{x_{i}\\}\to P$ of $k$-rational points approximating $P$,
if infinitely many points of $\\{x_{i}\\}$ are outside $B$ then
$\alpha(\\{x_{i}\\},P,L)\geq\gamma$.

  3. (c)

If $\alpha(P,L)<\gamma$ then $P\in B$ and $\alpha(P,L)=\alpha(P,L|_{B})$.

  4. (d)

If $P\in B$ and $\alpha(P,L|_{B})\geq\gamma$ then $\alpha(P,L)\geq\gamma$.

In other words, Theorem 2.4 says that if $L_{\gamma}$ is effective, then
$\alpha(P,L)\geq\gamma$, except maybe for sequences that lie in the image of
the stable base locus of $L_{\gamma}$.

The idea is to use Theorem 2.4 to show that most points don’t approximate $P$
very well, and then Theorem 2.3 to show that the terrific curve $C$ we will
find actually does approximate $P$ well.

##  3\. Main Theorem

Let $k$ be a number field, and let $v$ be a place of $k$. Let
$X\subset{\mathbb{P}}^{n}$ be a smooth cubic hypersurface defined over $k$,
and let $L$ be the hyperplane class on $X$.

If the dimension of $X$ is one, then of course there cannot be any line on
$X$, $k$-rational or not, so the claim is vacuous. Nevertheless, it is
classical111See for example the second theorem on page 98 of [Se]. that
$\alpha(P,L)=\infty$ for every point $P$ of $X$, and the  (Coba conjecture).
is vacuously true for $X$.

If the dimension of $X$ is two, then the question of computing $\alpha(P,L)$
is treated in [MR2]. In [MR2], the Picard group of the surface is assumed to
be entirely defined over $k$, but the argument for computing $\alpha(P,L)$ for
the anticanonical class – that is, the class that embeds $X$ as a cubic
hypersurface – depends only on the existence of a $k$-rational conic through
$P$, which is equivalent to the existence of a $k$-rational line on $X$ …
unless $P$ is the intersection of two conjugate irrational lines, in which
case the argument breaks down. We remedy this failing with the following
theorem, which also relaxes the smoothness hypothesis on the hypersurface.

######  Theorem 3.1.

Let $X\subset{\mathbb{P}}^{3}$ be a cubic surface (irreducible, but possibly
singular) defined over $k$. Let $P\in X(k)$ be any $k$-rational smooth point,
and let $L$ be the hyperplane class. Assume that there is a line $T$ on $X$
defined over $k$, and let $S_{P}$ be the intersection of $X$ with the tangent
plane at $P$. Assume further that the singular locus of $X$ is supported on
$T$. Then:

$\alpha(P,L)=\left\\{\begin{array}[]{ll}1&\mbox{if $P$ lies on a $k$-rational
line}\\\ 2&\mbox{if $P$ is an isolated point in $S_{P}(k)$}\\\ &\mbox{or if
$S_{P}$ has $k$-rational tangent lines at $P$}\\\
3/2&\mbox{otherwise}\end{array}\right.$

In all cases, there is a curve of best approximation to $P$.

Proof: Let $\pi\colon Y\to X$ be the blowup of $X$ at $P$, with exceptional
divisor $E$. Since $\pi^{*}L-2E$ is effective with stable base locus $S_{P}$,
Theorem 3.3 of [MR2] immediately shows that $\alpha_{L}(P)\geq 2$, except
possibly for sequences contained in $S_{P}$.

Now, $S_{P}$ is a plane curve of degree $3$ with a singularity at $P$, so we
can compute $\alpha(P,L|_{S_{P}})$ using Theorem 2.3. If $P$ is $v$-adically
isolated on $S_{P}$, then of course $\alpha(P,L|_{S_{P}})=\infty$, and so
$S_{P}$ is no threat to the computation of $\alpha$.

Otherwise, we can assume that $P$ is the $v$-adic limit of $k$-rational points
on $S_{P}$. Theorem 2.3 confidently lists the possibilities here:

  * •

If $S_{P}$ is reducible or non-reduced, then it’s the union of two
$k$-rational lines and $\alpha(P,L|_{S_{P}})=1$.

  * •

If $S_{P}$ is irreducible and there are two tangent lines to $S_{P}$ at $P$,
both defined over $k$, then $\alpha(P,L|_{S_{P}})=3$.

  * •

If $S_{P}$ is irreducible and there are two tangent lines to $S_{P}$ at $P$,
neither defined over $k$, then $\alpha(P,L|_{S_{P}})=3/2$. (Remember that $P$
is not $v$-adically isolated on $S_{P}$, so the tangent lines are defined over
the completion $k_{v}$.)

  * •

If $S_{P}$ is irreducible and there is only one tangent line to $S_{P}$ at
$P$, then $\alpha(P,L|_{S_{P}})=3/2$.

If you match up this list with the statement of the theorem, you hit paydirt:
the values of $\alpha(P,L|_{S_{P}})$ must equal the values of $\alpha(P,L)$
when they’re less than $2$, and we already knew $\alpha(P,L)\geq 2$ in all the
other cases.

It therefore remains only to show $\alpha(P,L)\leq 2$. In other words, we need
to find a curve $C$ such that $\alpha(P,L|_{C})\leq 2$.

It’s tempting to think “oh, I’ll just use the conic I get when I intersect the
plane through $P$ and $T$ with $X$”. That pretty useful, for sure, but what if
that “conic” is actually the union of two irrational lines? Problems then,
because that “conic” doesn’t actually give you an upper bound for
$\alpha(P,L)$.

So we roll up our sleeves and get to work. First off, if $X$ is singular along
all of $T$, then the hyperplane through $P$ and $T$ is the union of $T$ with a
line, so $\alpha(P,L|_{S_{P}})=1$. Thus, in what follows, we will assume that
$X$ has only isolated singularities, all supported on $T$ by assumption.

Let $f\colon\tilde{X}\to X$ be a minimal desingularization of $X$, let
$-K=f^{*}{\mathcal{O}}(1)$ be the anticanonical class on $\tilde{X}$, and let
$\tilde{T}$ be the class on $\tilde{X}$ of the strict transform of the
$k$-rational line $T$. It’s classical (see for example [BW]) that the
singularities of $X$ are canonical, so $f$ is a crepant resolution of $X$:
$K_{\tilde{X}}=f^{*}K_{X}$. This means that since $T.(-K_{X})=1$, we get
$\tilde{T}.(-K_{\tilde{X}})=1$ (curves contracted by $f$ intersect
$K_{\tilde{X}}$ trivially), and so by adjunction we have $\tilde{T}^{2}=-1$.

Let $A$ be the class $\tilde{T}-K$, so that $\deg_{-K}(A)=A.(-K)=4$ and
$A^{2}=4$. It is a straightforward calculation that $A$ is big and nef and
that $H^{0}(A)=5$. There is therefore a pencil
$\mathcal{F}\cong{\mathbb{P}}^{1}$ of curves in $|A|$ such that each curve in
$\mathcal{F}$ has a singularity at $P$. Moreover, adjunction tells us that a
smooth curve in $|A|$ has genus 1, so a general curve in $\mathcal{F}$ has
geometric genus zero with a double point at $P$. (Bertini’s Theorem tells us
that all but finitely many elements of $\mathcal{F}$ are irreducible.)

No two curves in $\mathcal{F}$ can have the same tangent line at $P$ because
$A^{2}=4$, and so the morphism
$\Psi\colon\mathcal{F}\to\mbox{Sym}^{2}(\mathcal{H})$ is injective.

(Um. Forgot to tell you about $\Psi$, sorry. $\mathcal{H}$ is the projective
line parametrizing tangent directions at $P$ in the Zariski tangent space at
$P$. $\Psi$ is given by $\Psi(C)=\\{t_{1},t_{2}\\}$, where $t_{1}$ and $t_{2}$
represent the two tangent directions to $C$ at $P$.)

But $\mathcal{F}$ and $\mathcal{H}$ are projective curves! So “$\Psi$ is
injective” means “$\Psi$ is surjective” too. Thus, every tangent direction at
$P$ must occur as a tangent direction of some curve in $\mathcal{F}$.

We will find a curve in $\mathcal{F}$ with irrational tangent lines that are
defined over $k_{v}$. Theorem 2.3 will tell us that this curve $C$ will be the
one we’re looking for, with $\alpha(P,L|_{C})=\deg C/2=4/2=2$.

If $B\subset\mbox{Sym}^{2}(\mathcal{H})$ is the image of $\Psi$, then its
preimage $\tilde{B}\subset\mathcal{H}\times\mathcal{H}$ is in the class
$(1,1)$, and therefore induces an involution on $\mathcal{H}$ whose quotient
is a surjective double cover $q\colon\mathcal{H}\to\mathcal{F}$. That’s a
fancy way of saying “$q$ maps (a pair of tangent directions) to (the curve in
$\mathcal{F}$ that they come from)”.

All we need now is a $k$-rational element $C\in\mathcal{F}$ such that
$q^{-1}(C)$ is the union of two distinct irrational points that are defined
over the completion $k_{v}$. In fact, we will show that there are infinitely
many of these.

To do this, it is enough to prove the following lemma, which to be honest is
probably not new:

######  Lemma 3.2.

Let $f\colon{\mathbb{P}}^{1}\to{\mathbb{P}}^{1}$ be a morphism of degree two,
defined over a number field $k$. For any place $v$ of $k$, there are
infinitely many points $P\in{\mathbb{P}}^{1}$ such that $f^{-1}(P)$ consists
of two points that are defined over $k_{v}$ but not $k$.

Proof of lemma: Let $Q\in{\mathbb{P}}^{1}(k)$ be a rational point where $f$ is
not ramified. Then there is an open $v$-adic neighbourhood $U$ of $Q$ such
that for all $R\in U$, $f$ is unramified at $R$.

The set $V=f(U)$ is therefore an open $v$-adic neighbourhood of the
$k$-rational point $f(Q)$ of ${\mathbb{P}}^{1}$, with the property that for
every point $P\in V$, the set $f^{-1}(P)$ consists of two points defined over
$k_{v}$. The set $V$ contains a positive proportion of the $k$-rational points
of ${\mathbb{P}}^{1}$ by height, but the Hilbert Irreducibility Theorem (see
for example section 9 of [Se]) implies that the proportion of $k$-rational
points of ${\mathbb{P}}^{1}$ by height whose preimage is two $k$-rational
points is zero. Therefore there are infinitely many points of $V$ whose
preimage is defined over $k_{v}$ but not $v$. $\clubsuit$

We therefore have a curve $C$ that is defined over $k$, irreducible, genus
zero, singular at $P$, and has tangent directions at $P$ that are not defined
over $k$, but are defined over the completion $k_{v}$. (Indeed, we have
infinitely many such curves.) By Theorem 2.6 of [MR2], $C$ contains a sequence
with approximation constant $\deg(C)/2=2$, and is therefore the curve we seek.

Thus, there is always a curve of best approximation to $P$, as desired, and
the value of $\alpha$ is as advertised. $\clubsuit$

Having treated the curve and surface cases, we will henceforth assume that the
dimension of $X$ is at least three.

######  Theorem 3.3.

Let $n\geq 4$ be an integer, $X\subset{\mathbb{P}}^{n}$ a smooth cubic
hypersurface defined over $k$, $P\in X(k)$ any $k$-rational point, and let $L$
be the hyperplane class. Assume that there is a line on $X$ defined over $k$,
and let $S_{P}$ be the intersection of $X$ with the tangent hyperplane at $P$.
Then:

$\alpha_{L}(P)=\left\\{\begin{array}[]{ll}1&\mbox{if $P$ lies on a
$k$-rational line}\\\ 2&\mbox{if $P$ is an isolated point in $S_{P}(k)$}\\\
3/2&\mbox{otherwise}\end{array}\right.$

In all cases, there is a curve of best approximation to $P$.

######  Remark 3.4.

Note that the three cases are mutually exclusive, as any $k$-rational line
through $P$ must lie on $S_{P}$.

Proof: First, notice that by Theorem 2.3, if $P$ lies on a $k$-rational line,
then $\alpha_{L}(P)=1$, and the line is a curve of best approximation. Thus,
we assume for the rest of the proof that $P$ does not lie on a $k$-rational
line.

Let $\pi\colon Y\to X$ be the blowup of $X$ at $P$, with exceptional divisor
$E$. Then $\pi^{*}L-2E$ is effective (the tangent hyperplane has multiplicity
2 at $P$), with stable base locus contained in $S_{P}$. By Theorem 2.4, this
means that $\alpha(P,L)\geq 2$, unless there is a sequence contained in
$S_{P}$ whose approximation constant is less than $2$.

We first consider the case where $P$ is an isolated rational point on $S_{P}$.
We just showed that $\alpha(P,L)\geq 2$. All that’s left is to show that
$\alpha(P,L)\leq 2$.

Let $H$ be any three-dimensional linear space that intersects the tangent
space at $P$ properly and transversely, contains $P$ and the line $\ell$, and
is defined over $k$. Let $V=H\cap X$ be the intersection of $H$ and $X$. By
Bertini’s Theorem we may choose $H$ so that $V$ is a cubic surface in $H$ that
is defined over $k$ and smooth at $P$, and whose singularities are supported
on $\ell$. Thus, by Theorem 3.1, we deduce that
$\alpha(P,L)\leq\alpha(P,L|_{V})\leq 2$.

We may therefore assume that $P$ is not an isolated rational point on $S_{P}$,
and restrict our attention to $S_{P}$. We will show that
$\alpha(P,L|_{S_{P}})=3/2$, and find a curve of best approximation to $P$ with
respect to $L|_{S_{P}}$. (Remember we assumed that $P$ doesn’t lie on a
$k$-rational line!)

Let $\phi\colon S_{P}\dashrightarrow{\mathbb{P}}^{n-2}$ be the linear
projection from the point $P$, considered inside the tangent hyperplane
$T_{P}\cong{\mathbb{P}}^{n-1}$. The fibres of $\phi$ (over ${\bar{k}}$) are
either single points or lines.

If $\phi$ is not dominant, therefore, then all the fibres are lines and
$S_{P}$ is a cone. Since $P$ is not an isolated rational point on $S_{P}$,
there is another rational point on $S_{P}$, and therefore a $k$-rational line
on $P$, which contradicts our earlier assumption.

Thus, we may assume that $\phi$ is dominant. In this case, we may change
coordinates in $T_{P}$ (which does not change $\alpha$) so that
$P=[0:0:\ldots:1]$ and the defining equation for $S_{P}$ is

$f(x_{1},\ldots,x_{n-1})+x_{n}g(x_{1},\ldots,x_{n-1})=0$

where $f$ is a homogeneous cubic and $g$ is a nonzero homogeneous quadric. The
rational map

$\psi({\mathbf{x}})=[x_{1}g(x_{1},\ldots,x_{n-1}):\ldots:x_{n-1}g(x_{1},\ldots,x_{n-1}):-f(x_{1},\ldots,x_{n-1})]$

is inverse to $\phi$, defined everywhere except $f=g=0$. Thus, $S_{P}$ is the
blowup of ${\mathbb{P}}^{n-2}$ along the scheme $f=g=0$, with the strict
transform of $g=0$ contracted.

More precisely, let $Z\subset T_{P}\times{\mathbb{P}}^{n-2}$ is the closure of
the graph of $\phi$. If $\pi_{i}$ is the projection of $Z$ onto the $i$th
factor, then

  * •

$\pi_{2}\colon Z\to{\mathbb{P}}^{n-2}$ is the blowup along $f=g=0$.

  * •

$\pi_{1}\colon Z\to S_{P}$ is an isomorphism away from the strict transform
$C$ of $g=0$ from $\pi_{2}$ and contracts $C$ to the point $P$.

Now let $\\{y_{i}\\}\subset S_{P}$ be any sequence of rational points
converging to $P$. We will show that $\alpha(\\{y_{i}\\},P,L)\geq 3/2$. The
idea is to relate the line bundle $L$ to the line bundle ${\mathcal{O}}(1)$ on
${\mathbb{P}}^{n-2}$, because we know all about $\alpha$ on projective space
already.

First, note that all of the $y_{i}$ lie in the subset of $S_{P}$ on which
$\phi$ is an isomorphism, since none of them is $P$, and none of the lines
through $P$ are $k$-rational. That means there is a unique sequence
$\\{z_{i}\\}$ of rational points on $Z$ such that $\pi_{1}(z_{i})=y_{i}\in
S_{P}$.

Now let’s pull back both $L$ and $H={\mathcal{O}}_{{\mathbb{P}}^{n-2}}(1)$ to
$Z$ and compare them there.

Let $L^{\prime}=L|_{S_{P}}$, a divisor class on $S_{P}$, and let
$B=\pi_{1}^{*}L^{\prime}$, a divisor class on $Z$.

Let $H={\mathcal{O}}_{{\mathbb{P}}^{n-2}}(1)$, a divisor class on
${\mathbb{P}}^{n-2}$.

Then $B=\pi_{2}^{*}H+C$ and $C=2\pi_{2}^{*}H-D$ for some effective divisor
$D$. (In case you forgot: $C$ is the strict transform of $g=0$.) Therefore
$2B-3C=2\pi_{2}^{*}H-C=D$, which is effective.

Therefore, up to multiplication by bounded functions compactly supported away
from zero, we have for any $\gamma<3/2$:

$\displaystyle H_{L^{\prime}}(y_{i})\operatorname{dist}(P,y_{i})^{\gamma}$
$\displaystyle=H_{B}(z_{i})\operatorname{dist}(P,y_{i})^{\gamma}$
$\displaystyle\geq H_{C}(z_{i})^{3/2}\operatorname{dist}(P,x_{i})^{\gamma}$
$\displaystyle\geq\operatorname{dist}(P,y_{i})^{-3/2}\operatorname{dist}(P,y_{i})^{\gamma}$
$\displaystyle\rightarrow\infty$

as $i\to\infty$. We conclude that $\alpha_{L^{\prime}}(\\{y_{i}\\},P)\geq
3/2$, and therefore $\alpha_{L^{\prime}}(P)\geq 3/2$, as desired.

All that’s left is to show that $\alpha_{L^{\prime}}(P)\leq 3/2$, which we can
do by exhibiting an approximating sequence that achieves $\alpha=3/2$. We can
do that by finding a curve through $P$ of $L^{\prime}$-degree 3 and a suitable
singularity at $P$.

Since $P$ is not an isolated rational point on $S_{P}$, it follows that the
quadric $g=0$ must have points defined over the completion $k_{v}$ of $k$ at
the place $v$. Let $Q$ be an irrational point of $g=0$ (which we may choose to
be away from the locus $f=0$), defined over a quadratic extension of $k$ that
is contained in $k_{v}$, and let $T$ be the line joining $Q$ to its conjugate.

The line $T$ is defined over $k$, so let $T^{\prime}$ be the closure of the
curve $\psi(T)$ in $S_{P}$. The map $\psi$ is an isomorphism away from $g=0$,
which is contracts to $P$, so $T^{\prime}$ is just $T$ with the point $Q$
glued to its conjugate. The $L^{\prime}$ degree of $T^{\prime}$ is 3, so
$T^{\prime}$ is a cubic curve $T^{\prime}$ with a node at $P$, with tangent
lines corresponding to $Q$ and its conjugate. Since the tangent lines to
$T^{\prime}$ at $P$ are defined over $k_{v}$ but not over $k$, by Theorem 2.3,
there is a sequence on $T^{\prime}$ with $\alpha=3/2$, and $T^{\prime}$ is a
curve of best approximation to $P$. $\clubsuit$

## References

  * [BD] Brandes, J. and Dietmann, R.; “Rational lines on cubic hypersurfaces”, Math. Proc. Cambridge Philos. Soc. 171 (2021), no. 1, 99–112.
  * [BW] Bruce, J.W.; Wall, C.T.C, “On the classification of cubic surfaces”, J. London Math. Soc. (2) 19 (1979), no. 2, 245–256.
  * [DW] Dietmann, R. and Wooley, T.; “Pairs of cubic forms in many variables”, Acta Arith. 110 (2003), no. 2, 125–-140.
  * [McK] McKinnon, David, “A conjecture on rational approximations to rational points”, J. Algebraic Geom., 16 (2007), 257–303.
  * [MR] McKinnon, D. and Roth, M., “Seshadri constants, Diophantine approximation, and Roth’s theorem for arbitrary varieties”, Invent. Math. (200), 513–583 (2015). DOI 10.1007/s00222-014-0540-1.
  * [MR2] McKinnon, D. and Roth, M., “An analogue of Liouville’s Theorem and an application to cubic surfaces”, Eur. J. Math. 2 (2016), no. 4, 929–959. DOI:10.1007/s40879-016-0113-5
  * [Se] Serre, J.-P., Lectures on the Mordell-Weil Theorem, Vieweg, 1997.


# PAT image reconstruction using augmented sparsity regularization with
practicable tuning of regularization weight

Nadaparambil Aravindakshan Rejesh
Imaging Systems Lab
Department of Electrical Engineering
Indian Institute of Science (IISc)
Bangalore 560012, India
<EMAIL_ADDRESS>
&Muthuvel Arigovindan
Imaging Systems Lab
Department of Electrical Engineering
Indian Institute of Science (IISc)
Bangalore 560012, India
<EMAIL_ADDRESS>
Corresponding author

###### Abstract

Among all tissue imaging modalities, photo-acoustic tomography (PAT) has been
getting increasing attention in the recent past due to the fact that it has
high contrast, high penetrability, and has capability of retrieving high
resolution. The reconstruction methods used in PAT plays a crucial role in the
applicability of PAT, and PAT finds particularly a wider applicability if a
model-based regularized reconstruction method is used. This is because such
methods can yield high quality images from measurements obtained using fewer
number of transducers. A crucial factor that determines the quality of
reconstruction in such methods is the choice of regularization weight.
Unfortunately, an appropriately tuned value of regularization weight varies
significantly with variation in the noise level, as well as, with the
variation in the high resolution contents of the image, in a way that has not
been well understood. There has been attempts to determine optimum
regularization weight from the measured data in the context of using
elementary and general purpose regularizations. In this paper, we develop a
practicable method for semi-automated tuning of the regularization weight in
the context of using a modern type of regularization that is specifically
designed for PAT image reconstruction. As a first step, we introduce a
relative smoothness constraint with a parameter; this parameter
computationally maps into the actual regularization weight, but, its tuning
does not vary significantly with variation in the noise level, and with the
variation in the high resolution contents of the image. Next, we construct an
algorithm that integrates the task of determining this mapping along with
obtaining the reconstruction. Finally we demonstrate experimentally that we
can run this algorithm with a nominal value of the relative smoothness
parameter—a value independent of the noise level and the structure of the
underlying image—to obtain good quality reconstructions.

##  1 Introduction

Photoacoustic tomography (PAT) has received an increasing interest over the
past two decades in soft tissue imaging since it combines the advantages of
ultrasound imaging and high contrast optical excitation [1, 2, 3, 4, 5, 6, 7].
This hybrid imaging modality has found various applications in biomedicine
such as brain imaging in small animals [8, 9], rheumatoid arthritis research
[10, 11] and breast imaging in humans [12, 13, 14]. In PAT, some of the
absorbed optical energy is converted to heat when the tissue is illuminated
with nanosecond laser pulses. The resulting thermoelastic expansion inside the
tissue creates an initial pressure rise; this propagates as ultrasound waves
towards the boundary where they are detected by an ultrasound transducer
array. These detected signals undergo a reconstruction step to get an estimate
of 2D image of initial pressure-rise, which represents the spatial
distribution of the physiological quantity under study. By using this
combination of optical absorption and acoustic wave propagation, PAT can
retrieve higher resolution with higher contrast from larger depths compared to
purely optical modalities. In most applications of PAT, the high contrast
stems from the contrast in the optical excitability of tissue constituents
such as haemoglobin. Hence the visualisation of vascular structures by
detecting signals from haemoglobin and oxyhaemoglobin becomes possible, which
can aid in more robust detection or diagnosis of disorders that are
characterised by abnormal proliferation of blood vessels.

The most common reconstruction methods employed are filtered backprojection
(FBP) [15, 16, 17, 18, 19] and time reversal [20, 21, 22, 23]. These standard
methods suffer from curved line artefacts and blurring, especially when the
noise level is high and/or when the placement of detectors for spatial
sampling of the detected photoacoustic (PA) signal is restricted [24]. In the
past decade, there has been of greater interest in model-based reconstruction
methods that solve the PAT reconstruction problem by using regularization. In
regularized method, the reconstruction is obtained as a minimizer of cost that
is a sum of a data-fitting term and a weighted regularization functional or
simply a regularization. The data fitting term measures the goodness of fit of
the candidate image to the measured data via a discrete forward model, and the
regularization imposes some prior belief on the required solution. The weight
applied on the regularization is called the regularization weight. Such
methods typically yield superior reconstructions even when the number of
transducer locations is limited [25, 26, 27, 28]. Depending on the prior
assumptions on the spatial characteristics of the image to be recovered, the
regularizers used for these reconstructions have many different forms [25, 26,
27, 28, 29, 30, 31, 32, 33, 34, 35, 36]. The traditional choice among these
regularization functionals is the Tikhonov regularization functional. This
functional is constructed as the sum of squares of image intensity or image
derivatives. Minimization of this leads to suppression of noise, and has been
effectively used for PAT image recovery [28, 29, 30, 37, 38, 39, 40, 41].
However minimization of squared derivatives leads to over-smoothing; to
alleviate this problem, several improvements were proposed [39, 40].
Exponential filtering of singular values in the framework of Tikhonov
filtering was shown to provide better quantitative accuracy with reduced bias
[41].

Later, sparse regularization approaches gained attention due to their ability
to perform well in scenarios where noise is high and/or number transducer
locations is low [42, 27, 43, 34, 44, 45, 46, 31, 32, 33, 47, 48, 49]. Early
sparsity based methods used wavelet regularization as the prior [42, 44].
Total variation (TV) regularization gained prominence as the demonstration of
its superiority over wavelet regularization was facilitated by the
developments in convex optimization [27, 31, 32, 33, 45, 46, 47, 48, 50]. The
success of TV regularization stems from the fact that most objects have
sparsely distributed derivative values. Among the TV based methods, first
order TV (TV-1) [27, 50, 32, 31, 46, 47, 48] uses first order derivatives in
their formulation, and have been demonstrated to recover the images with
abrupt edges. However, when the noise is high and/or number of transducers is
low, TV-1 regularization gives reconstructions in the form of piece-wise
constants. To alleviate this, second-order TV (TV-2) regularization, which
uses second order derivatives, was applied in many imaging inverse problems,
and Boink et al. applied TV-2 for PAT image recovery [33]. TV-2 regularization
produces reconstruction with reduce artefacts in the presence of high level of
noise and/or high level of under-sampling, and it has become one of the most
widely used regularization for solving imaging inverse problems.

Further developments includes total generalised variation (TGV) regularization
[51] which helps to combine the best of TV-1 and TV-2. It has been applied for
PAT image reconstruction [52]. Although TGV combines the best of TV-1 and TV-2
and hence superior to TV-2 for general image recovery, a recent study reported
by Boink et al [33] indicates that TGV’s performance is comparable to that of
TV-2 in the case of PAT image recovery. In our previous work [53], we
constructed a regularization functional that combines the intensity and the
second-order derivatives to exploit the fact that high intensities and high
second-order derivatives are jointly sparse in PAT images. We proposed a novel
optimization method to solve the reconstruction problem using this
regularization, and demonstrated that this regularization significantly
outperforms TV regularization. In this paper, we will call this regularization
as the augmented sparsity regularization or simply the augmented
regularization.

It is important to note that accurately tuning the regularization weight is
crucial for good quality reconstruction in model based method, and, in a
practical scenario, the regularization weight has to be tuned by only using
the measured data. Methods that determine the regularization weight by only
using the measured data are sometimes called automated tuning methods;
because, without these methods, the user has to choose the weight by trial-
and-error approach where the reconstruction is performed for a series of
values of the weight with a visual feed back on the quality of the
reconstructed image, which is called the manual tuning. Clearly, manual tuning
is tedious as the appropriate value of the weight vary drastically with
variations in the noise level and the structural contents of the underlying
image.

With our focus on automated tuning of regularization weight, we note that
Prakash et al. have reported method for automated tuning of regularization
weight for PAT image recovery using Tikhonov regularization [40, 54, 55, 30].
As sparsity based regularization retrieves better resolution, our specific
focus will be on tuning of regularization weight for PAT image recovery in the
context of using such regularizations. With this, we note that most sparsity
based methods proposed for PAT reconstruction use the test models themselves
for determining the optimal value for regularization. To be more specific,
these methods simulate measurements from test image models, and adjust the
regularization weights such that the reconstruction obtained from the
simulated measurements best matches with the model that generated the
measurements. Although this way of tuning for regularization weight is
unrealistic since the underlying model image is used for tuning, the goal in
these works has been to demonstrate the effectiveness of the regularization.
We will call this way of tuning for regularization weight as the oracle
tuning.

There are few dedicated sparsity methods that deal with the problem of
determining the optimal regularization weight from the measured images for
elementary inverse problems such as denoising and deblurring; they are based
on Stein’s Unbiased Risk Estimator (SURE) principle. SURE principle gives a
computational expression for means squares error between the reconstructed
image and the underlying image that generated data without involving the
latter in the expression; SURE based method determine the required
regularization weight by minimizing this estimate. Although this sounds
appealing, the computational expression is very complex; it is complex because
it involves derivatives of the reconstruction as a function, where the
derivatives are taken with respect to measured data values. Further, SURE
based methods requires knowing the exact values of noise parameter(s). Early
SURE methods were based on wavelet regularization [56, 57, 58]. As the total
variation based methods give better results, the focus moved on to total
variation based methods [59, 60]. The work of Ramani et al. [59] provides
general computational framework for SURE based determination of regularization
weight for various problems including total variation based denoising. The
work of Xue et al [60] allows determining the regularization weight for
solving deblurring using total variation regularization. Adopting this
approach for PAT problem is impractical since this method is based on
expensive Monte Carlo simulations.

In this paper, we develop a semi-automatic method for determining the
regularization weight from measured data. We do this through the following
steps:

  * •

We propose a novel smoothness measure, that is a monotonically decreasing
function of the regularization weight, $\lambda$. While appropriate value of
$\lambda$ varies drastically with variations in the noise level and structural
properties of the underlying image, this smoothness measure lies within a
narrow range if the value of $\lambda$ is chosen appropriately. We propose to
use this smoothness measure as the user parameter instead of using
regularization weight $\lambda$ directly.

  * •

We propose a method that jointly performs the reconstruction and determines
the corresponding value of $\lambda$ for a user defined value of the relative
smoothness measure. We consider both TV-2 regularization and a modified
version of the augmented regularization proposed in our recent work [53].

  * •

We demonstrate experimentally that we can run this algorithm with a fixed
nominal value of the relative smoothness parameter to obtain good quality
reconstructions independent of the noise level and the structure of the
underlying image. We compare the structural similarity (SSIM) scores of
reconstruction obtained this way to that of reconstructions obtained with
oracle tuning of the regularization weight, and show that the SSIM scores are
comparable. This means that, in a practical point of view, our work solves the
problem of determining the required regularization weight from measured
images.

We state that our method is semi-automatic, because, the above-mentioned
nominal value for the relative smoothness was determined by trial and error
method. However, it should be emphasized that a specific value works for all
data sets independent of noise level and the underlying image. Hence our
algorithm has practical significance although it is semi-automatic.

##  2 Detailed overview of derivative based regularization methods for model
based PAT image reconstruction

The key relation necessary for the development of regularized reconstruction
in PAT is expressing the discrete image formation model. To this end, the
first step is expressing the forward model for PAT image reconstruction in
terms of the relation between the laser energy deposited in the imaging
specimen, $E(\textbf{r},t)$, and pressure field, $p(\textbf{r},t)$. Here both
quantities are expressed as functions of space and time, with $\mathbf{r}$
representing 2D spatial location and $t$ representing time. This relation is
given by [24]

$\Big{(}\frac{\partial^{2}}{\partial
t^{2}}-c^{2}_{0}\;\nabla^{2}\Big{)}\;p(\textbf{r},t)=\Gamma\;\frac{\partial}{\partial
t}E(\textbf{r},t),$ (1)

Here $\Gamma$ is a dimensionless parameter called the Gr$\ddot{u}$neisen
coefficient, which describes the conversion efficiency of heat to pressure and
$c_{0}$ is the speed of sound in a homogeneous medium. By recognizing that the
temporal duration of the laser pulse is shorter than the temporal resolution
of the ultrasound detectors in most of the photoacoustic imaging applications,
the PA source $E(\textbf{r},t)$ may be approximated by
$E(\textbf{r})\delta(t)$, where $E(\textbf{r})$ is the density of deposited
energy. Then the solution to the differential equation (1) can be written as
[24]

$p(\textbf{r},t)=\cfrac{\Gamma}{4\pi c_{0}}\;\;\frac{\partial}{\partial
t}\;\int\limits_{|\textbf{r}-\textbf{r}^{\prime}|=c_{0}t}\;\dfrac{E(\textbf{r}^{\prime})}{|\textbf{r}-\textbf{r}^{\prime}|}\;d\textbf{r}^{\prime}.$
(2)

Here the initial pressure field ${p}_{0}(\textbf{r})=p(\textbf{r},t=0)$ can be
written as

${p}_{0}(\textbf{r})=\Gamma\;E(\textbf{r}).$ (3)

The initial pressure field represents the physical quantity of interest and
goal of PAT image reconstruction task is to recover ${p}_{0}(\textbf{r})$ from
the samples of $p(\textbf{r},t)$. The key point that makes this task
challenging is that ${p}_{0}(\textbf{r})$ is required to be determined for the
entire cross-sectional plane, whereas the samples of $p(\textbf{r},t)$ are
available only from the points lying in the periphery of the imaging specimen
where the transducers are located.

To perform the regularized reconstruction in a computer, one needs to have
discrete representation of the relation given above. To this end, the required
image itself is represented as a discrete 2D array, and the samples of
$p(\textbf{r},t)$ are represented in terms of numerical approximation of the
integral given above. Let $\\{\mathbf{r}_{i},i=0,\ldots,L-1\\}$ be the
transducer locations and let $\\{t_{j},j=0,\ldots,M_{t}-1\\}$ be the sampling
instants. Further, suppose ${p}_{0}(\textbf{r})$ is discretized in a
$N_{x}\times N_{y}$ cartesian grid, and suppose ${\bf p}_{0}$ represents
discrete samples of ${p}_{0}(\textbf{r})$ scanned into vector form with length
$N=N_{x}N_{y}$. Then the forward model can be represented by

$\textbf{m}=\textbf{H}\textbf{p}_{0}$ (4)

where ${\bf m}$ is the $LM_{t}\times 1$ vector with $\ell$th element
satisfying $\\{{\bf m}\\}_{\ell}=p({\bf r}_{i},t_{j})$ such that
$\ell=M_{t}i+j$, and ${\bf H}$ is the matrix representing the discrete
equivalent of the integral in the equation (2). Several methods to improve the
accuracy of the model by incorporating transducer responses and
heterogeneities in the medium have been proposed [31, 61, 62]. Including the
measurement noise, the modified imaging equation can be written as

$\textbf{m}=\textbf{H}\textbf{p}_{0}+{\boldsymbol{\eta}},$ (5)

where ${\boldsymbol{\eta}}$ represents measurement noise, which is Gaussian.

A reconstruction method in PAT aims to recover the initial pressure
distribution $\textbf{p}_{0}$ from the noisy transducer measurement data m.
When the number of transducers is less, the PA image reconstruction problem is
ill-posed, and hence constraints are imposed on the required solution in the
form of regularization. In this case, the image reconstruction problem can be
treated as an optimization problem where the solution is obtained by
minimizing a cost function. The reconstruction problem can be written as

$\hat{\textbf{p}}_{0}=\operatorname*{arg\,min}_{\textbf{x}}\;J(\textbf{x})$
(6)

where $J(\textbf{x})$ is the cost function and is given by

$J(\textbf{x})=||\textbf{m}-\textbf{H}\textbf{x}||_{2}^{2}+\lambda\;{R}(\textbf{x}).$
(7)

Here ${R}(\textbf{x})$ is the regularization functional and $\|\cdot\|_{2}$
represents the $l_{2}$ norm. The regularization parameter $\lambda$ controls
the amount of smoothness in the solution and fidelity to the measured data.
The regularization functional typically should be able to smooth the noise and
hence constructed using derivative terms.

Tikhonov regularization has been used in limited data cases [24, 48] and is
given by

${R}(\textbf{x})=\sum_{i}\|\textbf{D}_{o,i}\textbf{x}\|^{2}_{2}=\sum_{r=1}^{N}\sum_{i}\big{(}(\textbf{D}_{o,i}\textbf{x})_{r}\big{)}^{2}$
(8)

where $(\cdot)_{r}$ denotes the $r$th component of its vector argument, and
$\textbf{D}_{o,i}$ represents the matrix of $i^{th}$ derivative filter of
order $o$. For example, $\textbf{D}_{1,i},i=1,2$, are the matrix equivalents
of filtering by discrete filters that implement the operators
$\frac{\partial}{\partial x}$ and $\frac{\partial}{\partial y}$. Further,
$\textbf{D}_{2,i},i=1,2,3$, are the matrix equivalents of filtering by
discrete filters that implement the operators $\frac{\partial^{2}}{\partial
x^{2}},\frac{\partial^{2}}{\partial y^{2}}$ and
$\sqrt{2}\frac{\partial}{\partial x\partial y}$. The resulting minimization of
the convex quadratic cost function yields a closed form solution given by

$\hat{\textbf{p}}_{0}=\big{[}\textbf{H}^{t}\textbf{H}+\lambda\;\sum_{i}\textbf{D}_{o,i}^{t}\textbf{D}_{o,i}\big{]}^{-1}\textbf{H}^{t}\textbf{m}.$
(9)

As emphasized before, Tikhonov regularization precludes any large derivative
values and hence it tends to smooth edges in the reconstructed image. Note
that we use the matrix formulation for derivatives for gaining notational
convenience in describing the minimization method. The matrices,
$\textbf{D}_{o,i}$ represent filters that implement the derivatives
$\frac{\partial}{\partial x}$, $\frac{\partial}{\partial y}$,
$\frac{\partial^{2}}{\partial x^{2}}$, $\frac{\partial^{2}}{\partial y^{2}}$
and $\sqrt{2}\frac{\partial}{\partial x\partial y}$; these can be directly
applied to the images without building the matrices. On the other hand, the
data fidelity term requires building the matrix ${\bf H}$ explicitly. Total
Variation (TV) regularization [63, 46, 64, 31] can be expressed as

${R}_{TV}(\textbf{x})=\sum_{r=1}^{N}\sqrt{\sum_{i}\big{(}\textbf{D}_{o,i}\textbf{x}\big{)}_{r}^{2}}$
(10)

With $o=1$ we get TV-1 and with $o=2$ we get TV-2 regularization. The recent
development, TGV regularization can be expressed as

${R}_{TGV}(\textbf{x})=\underset{{\bf u}_{1},{\bf
u}_{2}}{\operatorname{argmin}}\;\;{\cal N}_{1}({\bf x},{\bf u}_{1},{\bf
u}_{2})+{\cal N}_{2}({\bf u}_{1},{\bf u}_{2})$ (11)

where

$\displaystyle{\cal N}_{1}({\bf x},{\bf u}_{1},{\bf u}_{2})=\sum_{r=1}^{N}$
$\displaystyle\left[((\textbf{D}_{1,1}\textbf{x}\big{)}_{r}-({\bf
u}_{1})_{r})^{2}\right.\left.+((\textbf{D}_{1,2}\textbf{x}\big{)}_{r}-({\bf
u}_{2})_{r})^{2}\right]^{0.5}$ (12) $\displaystyle{\cal N}_{2}({\bf
u}_{1},{\bf u}_{2})=\sum_{r=1}^{N}$
$\displaystyle\left[((\textbf{D}_{1,1}\textbf{u}_{1}\big{)}_{r})^{2}+((\textbf{D}_{1,2}\textbf{u}_{2}\big{)}_{r})^{2}+\right.\left.+0.5((\textbf{D}_{1,2}\textbf{u}_{1}\big{)}_{r}+(\textbf{D}_{1,1}\textbf{u}_{2}\big{)}_{r})^{2}\right]^{0.5}$
(13)

with ${\bf u}_{1}$ and ${\bf u}_{2}$ being $N\times 1$ vectors. Note that the
TGV regularization is expressed via minimization in $2N$ dimensional space.

The photoacoustic images have high contrast due to the differential absorption
of light in the near-infrared region by chromophores such as hemoglobin. Due
to this, high values of initial pressure, $p_{0}(\mathbf{r})$, are sparsely
distributed. Further, regions having high derivative values are also sparsely
distributed. This pattern was also observed for fluorescence images in the
work presented in [65], where the regularization was constructed by adding an
intensity term to second-order derivatives. The combined point-wise cost went
into a logarithmic function and summed over all pixels. Inspired by the
success of this form of regularization, we proposed a modification that suits
PAT reconstruction problem in [53]. We replaced the log by a fractional power
and expressed regularization as

${R}_{h}(\textbf{x},q)=\sum_{r=1}^{N}\left(\epsilon_{d}+\alpha\big{(}\textbf{x}\big{)}_{r}^{2}+(1-\alpha)\sum_{i}\big{(}\textbf{D}_{o,i}\textbf{x}\big{)}_{r}^{2}\right)^{q},$
(14)

where the weight $\alpha\in(0,1)$ controls the relative penalization, and
$\epsilon_{d}$ is a small positive real number. The advantage of this
modification is that it allowed an optimization strategy that can efficiently
handle non-convex cost function. We chose $q<0.5$ meaning that the resulting
cost functional is non-convex. We also proposed an optimization method for PAT
image reconstruction using the above regularization and demonstrated their
effectiveness using numerical experiments. We call ${R}_{h}(\textbf{x},q)$ as
the augmented sparsity regularization.

##  3 Proposed method

###  3.1 Cost function and reconstruction algorithm

Our goal is to make the augmented sparsity regularization more practicable by
constructing an iterative method that integrates the task of determining the
regularization weight and obtaining regularized reconstruction. This is
challenging, and to make this computationally tractable, we simplify the
augmented sparsity regularization by making it convex. Specifically, we
simplify the regularization given in the equation (14) as given below:

${R}_{hc}(\textbf{x},\alpha)=\sum_{r=1}^{N}\left(\alpha\big{(}\textbf{x}\big{)}_{r}^{2}+(1-\alpha)\sum_{i=1}^{3}\big{(}\textbf{D}_{2,i}\textbf{x}\big{)}_{r}^{2}\right)^{1/2}.$
(15)

Note that, we have also eliminated $\epsilon_{d}$ required for making the
functional differentiable, since non-differentiability can be easily handled
when the function is convex. With this, the required reconstruction is
expressed as

$\hat{\textbf{x}}=\operatorname*{arg\,min}_{\textbf{x}}\;J_{c}(\textbf{x},{\bf
H},{\bf m},\lambda,\alpha)$ (16)

where $J_{c}(\textbf{x},{\bf H},{\bf m},\lambda,\alpha)$ is the cost function
and is given by

$\begin{split}J_{c}(\textbf{x},{\bf H},{\bf
m},\lambda,\alpha)=&\frac{1}{n}\|\textbf{m}-\textbf{H}\textbf{x}\|_{2}^{2}+\lambda\;{R}_{hc}(\textbf{x},\alpha)+{\cal
B}_{[u]}({\bf x}),\end{split}$ (17)

with $n$ denoting the length of the measurement vector ${\bf m}$. Further,
${\cal B}_{[u]}(\cdot)$ is the cost function that enforces the candidate image
pixels ${\bf x}$ to be in the range $[0,u]$, where $u$ is the user-defined
upper bound. If any of the pixels of ${\bf x}$ is outside this range, ${\cal
B}_{[u]}(\cdot)$ becomes $\infty$, and it will have zero value otherwise. This
function imposes exact bound constraint on the required reconstruction. Note
that, in our original formulation given in [53], we only incorporated an
approximate bound constraint, and here we have the exact bound constraint.
This is again due to the fact that the exact bound constraint is easy to
implement since the overall cost function is convex. The above minimization
problem is very similar to the case where TV-2 functional is used instead of
${R}_{hc}(\textbf{x},\alpha)$. However, in order to develop the reconstruction
method that integrates the task of determining the regularization weight, we
need to write down a specific algorithm for minimizing $J_{c}(\cdot)$.

For minimizing this cost, we adopt alternating direction method of multipliers
(ADMM) approach [eck]. Implementation of ADMM requires separating the sub-
functionals of $J_{c}(\cdot)$ with auxiliary variables that are linked to the
primary variable ${\bf x}$ in terms of linear equations. To this end, we
define

${\bf D}_{\alpha}=\left[\begin{array}[]{c}\sqrt{\alpha}\;{\bf I}\\\
\sqrt{(1-\alpha)}\;{\bf D}_{2,1}\\\ \sqrt{(1-\alpha)}\;{\bf D}_{2,2}\\\
\sqrt{(1-\alpha)}\;{\bf D}_{2,3}\end{array}\right].$

With this, ${R}_{hc}(\textbf{x},\alpha)$ given in the equation can be
expressed as

${R}_{hc}(\textbf{x},\alpha)=\sum_{r=1}^{N}\left(\sum_{j=0}^{3}({\bf
D}_{\alpha}{\bf x})^{2}_{r+jN}\right)^{1/2}.$ (18)

To further simplify the notation, we introduce the interlaced slicing
operation given by

$({\bf v})_{[r]}=[({\bf v})_{r}\;({\bf v})_{r+N}\;({\bf v})_{r+2N}\;({\bf
v})_{r+3N}]^{t}.$

Note that ${\bf v}$ is in $\mathbb{R}^{4N}$ and $({\bf v})_{[r]}$ is in
$\mathbb{R}^{4}$ for $r=1,\ldots,N$. With this, we can express
${R}_{hc}(\textbf{x},\alpha)$ as

${R}_{hc}(\textbf{x},\alpha)=\sum_{r=1}^{N}\|({\bf D}_{\alpha}{\bf
x})_{[r]}\|_{2}$ (19)

For notational convenience in developing the ADMM algorithm, we define
following function for a vector ${\bf v}\in\mathbb{R}^{4N}$:

${\cal S}({\bf v})=\sum_{r=1}^{N}\|({\bf v})_{[r]}\|_{2}.$ (20)

With this, the minimization problem can be expressed as

$\displaystyle({\textbf{x}}^{*},{\bf b}^{*},{\bf d}^{*})=$
$\displaystyle\operatorname*{arg\,min}_{\textbf{x},{\bf d},{\bf
b}}\;\frac{1}{n}\|\textbf{m}-\textbf{H}\textbf{x}\|_{2}^{2}+\lambda{\cal
S}({\bf d})+{\cal B}_{[u]}({\bf b}),\;\;s.t.\;\;{\bf D}_{\alpha}{\bf x}={\bf
d},\;\;{\bf x}={\bf b}.$

The solution to the above constrained optimization problem can be obtained
through the so-called augmented Lagrangian which can be written as follows:

$\begin{split}L({\bf x},{\bf y},\hat{\bf
y})=&\frac{1}{n}\|\textbf{m}-\textbf{H}\textbf{x}\|_{2}^{2}+\bar{\cal S}({\bf
y})+\frac{\beta}{2}\|\bar{\bf D}_{\alpha}{\bf x}-{\bf
y}\|_{2}^{2}+\hat{y}^{t}\left(\bar{\bf D}_{\alpha}{\bf x}-{\bf
y}\right),\end{split}$ (21)

where $\bar{\bf D}_{\alpha}=\left[\begin{array}[]{c}{\bf I}\\\ {\bf
D}_{\alpha}\end{array}\right]$, and ${\bf y}=\left[\begin{array}[]{c}{\bf
b}\\\ {\bf d}\end{array}\right]$, and $\hat{\bf y}$ is Lagrange’s multiplier
vector with dimension $5N$. Here, $\beta$ is a penalty parameter, and
$\bar{\cal S}({\bf y})=\lambda{\cal S}({\bf d})+{\cal B}_{[u]}({\bf b})$.
Constrained optimization theory suggests that, with $\hat{\bf y}$ correctly
chosen, say $\hat{\bf y}=\hat{\bf y}^{*}$, minimizing $L({\bf x},{\bf
y},\hat{\bf y})$ with respect to ${\bf x}$ and ${\bf y}$ solves the
optimization problem given in the equation (3.1) [66]. This $\hat{\bf y}^{*}$
is the maximizer of $q(\hat{\bf y})$, which is the so-called the dual function
of the constrained optimization problem; this dual function is given by

$q(\hat{\bf y})=\min_{{\bf x},{\bf
y}}\frac{1}{n}\|\textbf{m}-\textbf{H}\textbf{x}\|_{2}^{2}+\bar{\cal S}({\bf
y})+\hat{y}^{t}\left(\bar{\bf D}_{\alpha}{\bf x}-{\bf y}\right).$ (22)

The method of ADMM allows to determine $\hat{\bf y}^{*}$ jointly with required
solution defined by

$({\textbf{x}}^{*},{\bf y}^{*})=\operatorname*{arg\,min}_{\textbf{x},{\bf
y}}L({\bf x},{\bf y},\hat{\bf y}^{*}),$

in terms of simple elementary minimizations without explicitly determining
$q(\hat{\bf y})$. Given a current estimate of the solution, say, $({\bf
x}^{(k)},{\bf y}^{(k)},\hat{\bf y}^{(k)})$, the next refined estimate in ADMM
is obtained as follows:

$\displaystyle\textbf{x}^{(k+1)}$
$\displaystyle=\operatorname*{arg\,min}_{\textbf{x}\in\mathbb{R}^{N}}L({\bf
x},{\bf y}^{(k)},\hat{\bf y}^{(k)}),$ (23) $\displaystyle\textbf{y}^{(k+1)}$
$\displaystyle=\operatorname*{arg\,min}_{\textbf{y}\in\mathbb{R}^{5N}}L({\bf
x}^{(k+1)},{\bf y},\hat{\bf y}^{(k)}),$ (24) $\displaystyle\hat{\bf
y}^{(k+1)}$ $\displaystyle=\hat{\bf y}^{(k)}+\beta(\bar{\bf D}_{\alpha}{\bf
x}^{(k+1)}-{\bf y}^{(k+1)}).$ (25)

The formula that gives the solution for first minimization problem can be
easily obtained because $L(\cdot,\cdot,\cdot)$ is quadratic with respect to
${\bf x}$. The equation that represents the minimum is given by

$\left({\bf H}^{t}{\bf H}+\beta\bar{\bf D}^{t}_{\alpha}\bar{\bf
D}_{\alpha}\right){\bf x}^{(k+1)}={\bf H}^{t}{\bf m}+\bar{\bf
D}^{t}_{\alpha}\left(\beta{\bf y}^{(k)}-\hat{\bf y}^{(k)}\right).$ (26)

The above equation can be solved in either of two ways. If the image to be
reconstructed is of moderate size, inverse of ${\bf H}^{t}{\bf
H}+\beta\bar{\bf D}^{t}_{\alpha}\bar{\bf D}_{\alpha}$ can be precomputed
outside the ADMM loop, and the required solution in the ADMM loop can be
computed as ${\bf x}^{(k+1)}=\left({\bf H}^{t}{\bf H}+\beta\bar{\bf
D}_{\alpha}^{t}\bar{\bf D}_{\alpha}\right)^{-1}\left({\bf H}^{t}{\bf
m}+\bar{\bf D}_{\alpha}^{t}\left(\beta{\bf y}^{(k)}-\hat{\bf
y}^{(k)}\right)\right)$. Otherwise, the equation (26) can be solved for ${\bf
x}^{(k+1)}$ using the method of conjugate gradients. Next, from the structure
of the function $L(\cdot,\cdot,\cdot)$, we deduce that the minimization
problem given in the equation (24) can be equivalently expressed as

$\textbf{y}^{(k+1)}=\operatorname*{arg\,min}_{\textbf{y}\in\mathbb{R}^{5N}}\;\;\bar{\cal
S}({\bf y})+\frac{\beta}{2}\|{\bf y}-\bar{\bf y}^{(k+1)}\|_{2}^{2}$ (27)

where $\bar{\bf y}^{(k+1)}=\bar{\bf D}_{\alpha}{\bf
x}^{(k+1)}+\frac{1}{\beta}\hat{\bf y}^{(k)}$. Note that $\bar{\cal S}({\bf
y})=\lambda{\cal S}({\bf d})+{\cal B}({\bf b})$ and ${\bf
y}=\left[\begin{array}[]{c}{\bf b}\\\ {\bf d}\end{array}\right]$. Let ${\bf
y}^{(k+1)}=\left[\begin{array}[]{c}{\bf b}^{(k+1)}\\\ {\bf
d}^{(k+1)}\end{array}\right]$ and $\bar{\bf
y}^{(k+1)}=\left[\begin{array}[]{c}\bar{\bf b}^{(k+1)}\\\ \bar{\bf
d}^{(k+1)}\end{array}\right]$ denote the corresponding partitions. Then the
minimization problem of equation (27) can be written as the following sub-
problems:

$\displaystyle\textbf{d}^{(k+1)}$
$\displaystyle=\operatorname*{arg\,min}_{\textbf{d}\in\mathbb{R}^{4N}}\;\;\lambda{\cal
S}({\bf d})+\frac{\beta}{2}\|{\bf d}-\bar{\bf d}^{(k+1)}\|_{2}^{2}.$ (28)
$\displaystyle\textbf{b}^{(k+1)}$
$\displaystyle=\operatorname*{arg\,min}_{\textbf{b}\in\mathbb{R}^{N}}\;\;{\cal
B}({\bf b})+\frac{\beta}{2}\|{\bf b}-\bar{\bf b}^{(k+1)}\|_{2}^{2}.$ (29)

The solution to the second sub-problem is well-known and can be represented by
$\textbf{b}^{(k+1)}={\cal P}(\bar{\bf b}^{(k+1)})$ where ${\cal P}(\cdot)$.
denotes the clipping of the components of the vector by the range $[0,u]$.
Next, we note that the first sub-problem is separable with respect to the
slicing operation $(\cdot)_{[r]}$. In other words, the first subproblem can be
written as

$(\textbf{d}^{(k+1)})_{[r]}=\operatorname*{arg\,min}_{\textbf{v}\in\mathbb{R}^{4}}\lambda\|{\bf
v}\|_{2}+\frac{\beta}{2}\left\|{\bf v}-(\bar{\bf
d}^{(k+1)})_{[r]}\right\|_{2}^{2}.$ (30)

This means that we can compute $(\textbf{d}^{(k+1)})_{[r]}$ for $r=1,\ldots,N$
and assemble to get the solution $\textbf{d}^{(k+1)}$. The solution to the
problem of equation (30) is well known and can be expressed as [67]
$(\textbf{d}^{(k+1)})_{[r]}=\textbf{Prox}_{l_{2}}((\bar{\bf
d}^{(k+1)})_{[r]},\lambda/\beta)$, where $\textbf{Prox}_{l_{2}}({\bf
z},\lambda)$ is given by

$\textbf{Prox}_{l_{2}}({\bf z},\lambda)=\frac{\max(\|{\bf
z}\|_{2}-\lambda,0)}{\|{\bf z}\|_{2}}{\bf z}$ (31)

Considering the choice of $\beta$, it is sufficient to set it equal to $1$ for
all measured datasets, although minor improvement in the convergence speed can
be obtained by adjusting the value.

###  3.2 Modified algorithm with relative smoothness constraint

To introduce the idea of relative smoothness constraint, let ${\bf H}_{f}$ be
the given measurement matrix and let ${\bf H}$ be the matrix obtained by
removing a fraction of rows equal to $\Delta$. In other words,
$n/n_{f}=1-\Delta$ where $n$ and $n_{f}$ are number of rows in the matrices
${\bf H}$ and ${\bf H}_{f}$ respectively. Let ${\bf m}_{f}$ and ${\bf m}$ be
the corresponding measurement vectors. Let ${\bf x}(\lambda)$ denote the
reconstruction obtained from ${\bf m}$ by using the regularization weight
$\lambda$. In other words, let

$\displaystyle{\textbf{x}}(\lambda)=$
$\displaystyle\operatorname*{arg\,min}_{\textbf{x}}\;\;J_{c}(\textbf{x},{\bf
H},{\bf m},\lambda,\alpha).$ (32)

where

$\displaystyle J_{c}(\textbf{x},{\bf H},{\bf m},\lambda,\alpha)=$
$\displaystyle\frac{1}{n}\|\textbf{m}-\textbf{H}\textbf{x}\|_{2}^{2}+\lambda\;{R}_{hc}(\textbf{x},\alpha)+{\cal
B}_{[u]}({\bf x}).$

Consider

$\displaystyle J_{c}(\textbf{x},{\bf H}_{f},{\bf m}_{f},\lambda,\alpha)=$
$\displaystyle\frac{1}{n_{f}}\|\textbf{m}_{f}-\textbf{H}_{f}\textbf{x}\|_{2}^{2}+\lambda\;{R}_{hc}(\textbf{x},\alpha)+{\cal
B}_{[u]}({\bf x}).$

Now, our hypothesis is the following: if $\lambda$ is chosen such that ${\bf
x}(\lambda)$ does not fit noise, then the data fidelity terms
$\frac{1}{n}\|\textbf{m}-\textbf{H}{\bf x}(\lambda)\|_{2}^{2}$ and
$\frac{1}{n_{f}}\|\textbf{m}_{f}-\textbf{H}_{f}{\bf x}(\lambda)\|_{2}^{2}$
will also be nearly identical because of the normalization with number of
measurements. However, in practice, the difference
$\left\|\frac{1}{n}\|\textbf{m}-\textbf{H}{\bf
x}(\lambda)\|_{2}^{2}-\frac{1}{n_{f}}\|\textbf{m}_{f}-\textbf{H}_{f}{\bf
x}(\lambda)\|_{2}^{2}\right\|$ will be monotonically decreasing. As a result,
the difference $|J_{c}({\bf x}(\lambda),{\bf H},{\bf
m},\lambda,\alpha)-J_{c}({\bf x}(\lambda),{\bf H}_{f},{\bf
m}_{f},\lambda,\alpha)|$ will also be monotonically decreasing w.r.t.
$\lambda$. Based on these observations, we define the following:

${S}(\lambda)=\frac{|{\cal J}_{f}({\textbf{x}}(\lambda),\lambda)-{\cal
J}({\textbf{x}}(\lambda),\lambda)|}{0.5({\cal
J}_{f}({\textbf{x}}(\lambda),\lambda)+{\cal
J}({\textbf{x}}(\lambda),\lambda))},$ (33)

where ${\cal J}_{f}({\textbf{x}}(\lambda),\lambda)=J_{c}({\bf x}(\lambda),{\bf
H}_{f},{\bf m}_{f},\lambda,\alpha)$ and ${\cal
J}({\textbf{x}}(\lambda),\lambda)=J_{c}({\bf x}(\lambda),{\bf H},{\bf
m},\lambda,\alpha)$. We propose that the user can specify the required
smoothing indrectly as an upper bound on ${S}(\lambda)$, say $\epsilon$, and
the required value for $\lambda$ can be determined as the lowest value of
$\lambda$ that satisfy ${S}(\lambda)\leq\epsilon$. The main advantage is that
${S}(\lambda)$ is a measure of adeqaucy of the magnitude of $\lambda$ for a
given measured data set, and variation in the reconstructed image quality with
respect to changes in ${S}(\lambda)$ will be much narrower than the variation
in the reconstructed image quality with respect to changes in $\lambda$.

Let ${\lambda}(\epsilon)$ be the lowest value of $\lambda$ that satisfy
${S}(\lambda)\leq\epsilon$. In other words, let
${\lambda}(\epsilon)=min\\{\lambda:{S}(\lambda)\leq\epsilon\\}$. Since
$S(\lambda)$ is monotonic, this is the same as the condition
${S}(\lambda(\epsilon))=\epsilon$. Now, since solving this exactly is
impractical, we use the following strategy: we define monotonically increasing
sequence $\\{\lambda_{i},i=0,1,2,\ldots\\}$ such that
$S(\lambda_{0})>\epsilon$, and find a index $m$ such that
$S(\lambda_{m-1})>\epsilon$ and $S(\lambda_{m})\leq\epsilon$. We then declare
$\lambda_{m}$ as an approximation for ${\lambda}(\epsilon)$. We define the
sequence in the form $\lambda_{i}=\lambda_{0}k^{i}$ where $k$ is the real
number of the form $k=1+\delta$ with $\delta$ being a small positive number.
Here, since $k$ determines how well $\lambda_{m}$ approximates
${\lambda}(\epsilon)$, we denotes this approximation by
${\lambda}_{k}(\epsilon)$. Determining ${\lambda}_{k}(\epsilon)$ this way and
obtaining the corresponding reconstruction can be represented by the schematic
given in Figure 1. In summary, the algorithm of Figure 1 performs
reconstruction without the user having to specify the regularization weight
directly, and it jointly determines the required regularization weight from
the user-specified relative smoothness $\epsilon$. One apparent problem in the
algorithm is that it brings two additional parameters, viz., $\Delta$ and $k$.
However, as we will demonstrate in the experiments section, these parameters
can be set at a nominal value independent of the measurement vector ${\bf m}$.
In the schematic of Figure 1, we have given computationally efficient formula
for ${S}(\lambda)$ by exploiting the fact that the terms
$J_{c}({\textbf{x}}(\lambda),{\bf H}_{f},{\bf m}_{f},\lambda,\alpha)$ and
$J_{c}({\bf x}(\lambda),{\bf H},{\bf m},\lambda,\alpha)$ differ only by the
data fitting term. We also used the fact that ${\cal B}_{[u]}(\cdot)$ will be
zero at ${\bf x}(\lambda)$. We denote the overall operation of the
reconstruction and tracking given in Figure 1 by $({\bf
x}^{*},\lambda_{k}(\epsilon))\leftarrow{\cal T}({\bf x}_{0},\lambda_{0},k,{\bf
H},{\bf H}_{f},{\bf m},{\bf m}_{f},\alpha)$. Note that, here, ${\bf x}_{0}$
denote the initial candidate for reconstruction and $\lambda_{0}$ is the
starting value for used for tacking $\lambda$.

Figure 1: Reconstruction with regularization weight tracking.
$\frac{n}{n_{f}}=1-\Delta$ where $n$ and $n_{f}$ are number of elements in
${\bf m}$ and ${\bf m}_{f}$ respectively.

The above scheme only represents an intermediate stage towards developing more
pratical method, which is our next goal. The main drawback in the above scheme
is that the tracking uses full minimization of the cost for each candidate
value of $\lambda$. In the remainder of this section, we develop a more
efficient approach that integrates the task of determining an approximation
for $\lambda_{k}(\epsilon)$ from a given $\epsilon$ and the task of obtaining
the reconstruction. To this end, we introduce an approximation for
${S}(\lambda)$, denoted by $\tilde{S}(\lambda)$, as given below:

$\tilde{S}(\lambda)=\frac{|\hat{\cal
J}_{f}(\tilde{\textbf{x}}(\lambda),\lambda)-\hat{\cal
J}(\tilde{\textbf{x}}(\lambda),\lambda)|}{0.5(\hat{\cal
J}_{f}(\tilde{\textbf{x}}(\lambda),\lambda)+\hat{\cal
J}(\tilde{\textbf{x}}(\lambda),\lambda))},$ (34)

where $\tilde{\textbf{x}}(\lambda)$ denotes a partial reconstruction obtained
by performing a few ADMM iterations, and $\hat{\cal J}_{f}(\cdot)$ and
$\hat{\cal J}(\cdot)$ denote the costs obtained from ${\cal J}_{f}(\cdot)$ and
${\cal J}(\cdot)$ by removing ${\cal B}_{[u]}(\cdot)$. The main assumption
that will lead the proposed algorithm is that $\tilde{S}(\lambda)$ will have a
similar type of behaviour provided that the number of ADMM iterations is
sufficiently large. Based on this assumption, we propose a modification to the
scheme of Figure 1, by merging the minimization of $J_{c}(\cdot)$ with
tracking for $\lambda$. As first step, we define the modified tracking denoted
by $(\hat{\bf x},\hat{\lambda}_{k}(\epsilon))\leftarrow{\cal T}_{(M)}({\bf
x}_{0},\lambda_{0},k,{\bf H},{\bf H}_{f},{\bf m},{\bf m}_{f},\alpha,\epsilon)$
that is obtained by replacing the minimization ${\bf
x}(\lambda)\leftarrow\underset{\bf x}{\operatorname{argmin}}\;J_{c}({\bf
x},{\bf H},{\bf m},\lambda,\alpha)$ by $M$ cycles of ADMM iterations, denoted
by $\tilde{\bf x}(\lambda)\leftarrow ADMM_{(M)}({\bf x},{\bf H},{\bf
m},\lambda,\alpha)$. Here ${\bf x}$ in $ADMM_{(M)}({\bf x},{\bf H},{\bf
m},\lambda,\alpha)$ denotes the starting point of the iterations. This scheme
is given in Figure 2. Note that this approximate scheme brings another
parameter, $M$, the number of ADMM iterations performed before each
computation of $\tilde{S}(\lambda)$. Again, as we will demonstrate, this can
also be fixed at a nominal value independent of the measurement vector ${\bf
m}$.

The final proposed algorithm is the repeated application of this tracking with
reconstruction, ${\cal T}_{(M)}(.)$, where the pair $(\hat{\bf
x},\hat{\lambda}_{k}(\epsilon))$ returned by each tracking becomes the
initialisation to the next tracking. This process is repeated until relative
change in the reconstructed image is lower than the user-defined tolerance
$\epsilon_{t}$. This is given in Figure 3. Comparing with the original
tracking scheme of Figure 1, there is no guarantee that this approximate
scheme will produce the same results. However, as we demonstrate in the
experiment section, the result produced by the approximate scheme has a
quality that is comparable to that of the reconstruction obtained by oracle
tuning. Note that the final parameter, $\epsilon_{t}$, is the termination
tolerance. As any iterative optimization based image reconstruction, this
parameter is set to small positive real number.

Figure 2: Regularization weight tracking with partial reconstruction Figure 3:
Overall proposed algorithm. ${\cal T}({\bf x},\lambda,k,{\bf H},{\bf
H}_{f},{\bf m},{\bf m}_{f},\alpha,\epsilon)$ represents the algorithm of
Figure 2

##  4 Experimental results

Figure 4: Numerical phantoms used for evaluation of the proposed method

Four numerical phantoms viz. BloodVessel, Derenzo, PAT, and RatBrain, given in
Figure 4 were used to evaluate PAT reconstruction methods. The first three are
synthetic phantoms, whereas the last was derived from the MRI angiogram of rat
brain published by Pastor et al [pastor2017general]. All were normalized to
the size 128 × 128 with a corresponding physical size of 12.8 mm × 12.8 mm as
models for generating the synthetic data. The measurement data was generated
using a numerical approximation of forward model described in [68] and added
with Gaussian noise to form simulated data having SNR levels of 15 dB, 20 dB,
25 dB and 30 dB. The number of transducers to mimic the limited data scenario
in our experiments was set to 16. The sound speed in the medium was assumed to
be 1.5 mm/$\mu$ s, and we considered the medium to be homogeneous with no
dispersion or absorption of sound. We consider three forms of reconstruction
methods for comparative evaluation:

  * •

AR-TR: Reconstruction using augmented regularization with proposed tracking of
$\lambda$, i.e., the method given in Figure 3.

  * •

AR-O: Reconstruction using augmented regularization with oracle tuning of
$\lambda$

  * •

TV2-O: Reconstruction using TV-2 regularization with oracle tuning of
$\lambda$

Our primary goal is to demonstrate the power of augmented regularization with
the regularization weight determined using the proposed tracking method of
Figure 3. This is taken care by the above list in two ways: (i) comparing AR-
TR with AR-O will show how well the power of augmented regularization will be
accessible in a practical scenario where the regularization weight is
determined from the measurement vector using the proposed tracking; it will
also quantify the loss incurred by using the approximations involved in the
proposed tacking; (ii) comparing AR-TR with TV2-O will demonstrate that
augmented regularization will still have its advantage over the best
performing regularization in the literature (TV-2) despite the above-mentioned
loss.

The bound on the relatives smoothness, $\epsilon$, was fixed at $0.06$
uniformly for all measurement sets independent of noise level and independent
of test model used for simulating measurements. Further, the other parameter
involved in the tracking algorithm were set as follows: $\Delta=0.1$,
$k=1.05$, and $M=50$. Again these settings were kept fixed for all input data
sets. The relative weight $\alpha$ was kept fixed at $0.5$ uniformly for all
data sets as we did in our original work [53]. Finally, the termination
tolerance $\epsilon_{t}$ was set to $10^{-4}$. The results are presented in
Table 1. We observe the following from the tables:

  * •

SSIM scores of reconstructions obtained from AR-TR are comparable to that of
the reconstructions obtained from AR-O; in most of cases, the difference is in
the third decimal point, and in rare cases, the difference is in the second
decimal point. This is significant because AR-O uses the model to determine
the value for $\lambda$ whereas AR-TR does not.

  * •

SSIM scores of reconstructions obtained from AR-TR are significantly better
than that of TV2-O, although AR-TR does not use the test model for determining
the value for $\lambda$ as opposed to TV2-O, which uses the model to determine
the value for $\lambda$.

By noting the fact that the required $\lambda$ values for AR-TR were tracked
from a fixed value of $\epsilon$, we conclude from the first observation that
the proposed method provide a practical solution to the problem of determining
the regularization weight for PAT image reconstruction. The second observation
imply that the advantage of augmented regularization over the standard TV
regularization is accessible via a practical method for tuning $\lambda$. Note
that in our previous work where we introduced the augmented regularization, we
demonstrated its superiority over the standard TV only by using oracle tuning.

Note that, among the published methods for regularization weight tuning with
total variation regularization, the method of Xue et al. [60] is closest
competitor for the the proposed method, AR-TR. However, comparison with TV2-O
eliminates the need to compare with method of Xue et al. because TV2-O is
bound to give reconstruction that is better than that of the method of Xue et
al. due to the oracle tuning in TV2-O. Besides, the method of Xue et al. is
based on Monte Carlo simulations, and implementing it for PAT reconstruction
will be impractical as PAT forward model is more complex than the blurring
model considered by Xue et al. In the remainder of this section, we will
display images corresponding to two selected cases from Table 1 and Table 2.

Phtm. | SNR (dB) | AR-TR | TV2-O | AR-O
---|---|---|---|---
RatBr | 15 | .916 | .791 | .919
20 | .955 | .877 | .957
25 | .977 | .925 | .979
30 | .986 | .956 | .988
BldVes | 15 | .919 | .719 | .928
20 | .952 | .827 | .964
25 | .980 | .881 | .983
30 | .993 | .945 | .993
Derz | 15 | .995 | .959 | .997
20 | .999 | .977 | .999
25 | .999 | .986 | .999
30 | .999 | .989 | .999
PAT | 15 | .985 | .946 | .985
20 | .990 | .962 | .990
25 | .990 | .973 | .992
30 | .991 | .978 | .994
Table 1:  SSIM scores of reconstruction with varying input noise levels with
16 transducers and relative smoothness parameter ($\epsilon$) value of 0.06.
Figure 5: Comparison of reconstructions obtained from simulated data with 16
transducers and 15 dB measurement noise. (A): RatBrain phantom model (the
maximum initial Pressure rise is assumed to be 1Pa); (B): Reconstruction
obtained by AR-TR; (C) Reconstruction obtained TV2-O; (D): Reconstruction
obtained by AR-O.

Figure 5 compares the reconstructed results for the RatBrain phantom with the
measurement set corresponding to SNR level of 15 dB. The original image is
given in Figure 5.A while Figure 5.B, 5.C and 5.D gives the reconstructed
images from AR-TR, TV2-O, and AR-O respectively. Clearly, AR-TR outperforms
TV2-O, which is also confirmed by the fact that the SSIM score of AR-TR’s
reconstruction is 0.125 higher than that of TV2-O’s reconstruction. Further,
AR-TR’s reconstruction has almost the same quality as that of AR-O. This is
confirmed by the fact that the SSIM score of AR-O’s reconstruction is only
0.003 higher than that of AR-TR’s reconstruction. Figure 6 shows the scan-line
based intensity profile of the reconstructions of Figure 5 . The intensity
profiles show that the proposed method, AR-TR, follows the ground truth as
good as the oracle tuning method, AR-O, which is also confirmed by the small
difference in the SSIM scores.

Figure 6: Scan line based intensity profiles of reconstructed images from
Fig:5.  Figure 7: Comparison of reconstructions obtained from simulated data
with 16 transducers and 15 dB measurement noise. (A):BloodVessel phantom model
(the maximum initial Pressure rise is assumed to be 1Pa); (B): Reconstruction
obtained by AR-TR; (C) Reconstruction obtained TV2-O; (D): Reconstruction
obtained by AR-O.

Next, for another closer view, we select the reconstruction corresponding to
measured data with 15 dB SNR obtained using BloodVessel phantom. The images
are given in Figure 7. Figure 7.A shows the BloodVessel phantom and Figures
7.B, 7.C and 7.D gives the reconstructions from the AR-TR, TV2-O, and AR-O
respectively. Here also, AR-TR is able to produce a reconstruction with a
quality that is comparable to the quality of the reconstructions produced by
AR-O, which is confirmed by a mere difference of 0.009 in SSIM. Further, AR-TR
retains the advantage that augmented regularization has over the TV-2
regularization, which is confirmed by the improvement of 0.2 in SSIM score.
Figure 8 shows the scan-line based intensity profile of the reconstructions in
Figure 7. Here also, our method closely follows the ground truth as good as
the oracle tuning method, which is also confirmed by the small difference in
the SSIM scores.

Figure 8: Scan line based intensity profiles of reconstructed images from
Fig:7.

##  5 Conclusion

Currently, there is no method for determining the required regularization
weight from measured data for PAT image reconstruction with sparsity
regularization; users of a sparsity method have to determine the weight by
trail and error approach where the reconstruction is performed for a series of
values for the regularization weight, $\lambda$, with a visual feed back on
the quality of the reconstructed image. This is impractical because the
appropriate value of $\lambda$ vary drastically with variations in the noise
level and the structural contents of the underlying image. In this paper, we
addressed this problem by introducing a novel relative smoothness parameter as
an implicit function of $\lambda$ that is defined through a partial
reconstruction. We constructed an iterative method that jointly determines,
for given value of relative smoothness, the corresponding $\lambda$ and the
required reconstruction with augmented sparsity regularization. We
demonstrated experimentally that we can get good quality reconstructions from
this new method by setting the relativeness parameter at a fixed nominal value
independent of the noise level and the nature of the underlying image. We
showed that the SSIM scores of reconstructions obtain this way were close to
the SSIM scores of reconstructions obtained from the scheme where the
underlying images themselves were used for determining the correct value for
$\lambda$. Further, other parameters involved in the definition of the
relative smoothness measured also have been set at fixed values independent of
noise level and the underlying image. This means that, in practical point of
view, our method solves the problem of determining the regularization weight
from measured data.

## References

  * [1] Wang L V and Yao J 2016 Nature methods 13 627
  * [2] Pramanik M, Ku G, Li C and Wang L V 2008 Medical physics 35 2218–2223
  * [3] Beard P 2011 Interface focus 1 602–631
  * [4] Zhou Y, Yao J and Wang L V 2016 Journal of biomedical optics 21 061007
  * [5] Upputuri P K and Pramanik M 2016 Journal of biomedical optics 22 041006
  * [6] Li L, Zhu L, Ma C, Lin L, Yao J, Wang L, Maslov K, Zhang R, Chen W, Shi J et al. 2017 Nature biomedical engineering 1 1–11
  * [7] Han Y, Ding L, Ben X L D, Razansky D, Prakash J and Ntziachristos V 2017 Optics letters 42 979–982
  * [8] Xia J and Wang L V 2013 IEEE Transactions on Biomedical Engineering 61 1380–1389
  * [9] Deán-Ben X, Gottschalk S, Mc Larney B, Shoham S and Razansky D 2017 Chemical Society Reviews 46 2158–2198
  * [10] van den Berg P J, Daoudi K, Moens H J B and Steenbergen W 2017 Photoacoustics 8 8–14
  * [11] Jo J, Tian C, Xu G, Sarazin J, Schiopu E, Gandikota G and Wang X 2018 Photoacoustics 12 82–89
  * [12] Manohar S and Dantuma M 2019 Photoacoustics 16 100134
  * [13] Toi M, Asao Y, Matsumoto Y, Sekiguchi H, Yoshikawa A, Takada M, Kataoka M, Endo T, Kawaguchi-Sakita N, Kawashima M et al. 2017 Scientific reports 7 1–11
  * [14] Schoustra S M, Piras D, Huijink R, Op’t Root T J, Alink L, Kobold W M F, Steenbergen W and Manohar S 2019 Journal of biomedical optics 24 121909
  * [15] Finch D and Patch S K 2004 SIAM journal on mathematical analysis 35 1213–1240
  * [16] Kunyansky L A 2007 Inverse problems 23 373
  * [17] Xu M and Wang L V 2002 IEEE transactions on medical imaging 21 814–822
  * [18] Xu M and Wang L V 2005 Physical Review E 71 016706
  * [19] Haltmeier M, Schuster T and Scherzer O 2005 Mathematical methods in the applied sciences 28 1919–1937
  * [20] Xu Y and Wang L V 2004 Physical review letters 92 033902
  * [21] Burgholzer P, Matt G J, Haltmeier M and Paltauf G 2007 Physical Review E 75 046706
  * [22] Hristova Y, Kuchment P and Nguyen L 2008 Inverse problems 24 055006
  * [23] Treeby B E and Cox B T 2010 Journal of biomedical optics 15 021314
  * [24] Rosenthal A, Ntziachristos V and Razansky D 2013 Current Medical Imaging 9 318–336
  * [25] Buehler A, Rosenthal A, Jetzfellner T, Dima A, Razansky D and Ntziachristos V 2011 Medical physics 38 1694–1704
  * [26] Qin T, Zheng Z, Zhang R, Wang C and Yu W 2019 Physics in Medicine & Biology 64 195004
  * [27] Yao L and Jiang H 2011 Biomedical optics express 2 2649–2654
  * [28] Sanny D R, Prakash J, Kalva S K, Pramanik M and Yalavarthy P K 2018 Journal of biomedical optics 23 100502
  * [29] Rejesh N A, Pullagurla H and Pramanik M 2013 JOSA A 30 1994–2001
  * [30] Shaw C B, Prakash J, Pramanik M and Yalavarthy P K 2013 Journal of Biomedical Optics 18 080501
  * [31] Huang C, Wang K, Nie L, Wang L V and Anastasio M A 2013 IEEE transactions on medical imaging 32 1097–1110
  * [32] Arridge S R, Betcke M M, Cox B T, Lucka F and Treeby B E 2016 Inverse Problems 32 115012
  * [33] Boink Y E, Lagerwerf M J, Steenbergen W, van Gils S A, Manohar S and Brune C 2018 Physics in Medicine & Biology 63 045018
  * [34] Frikel J and Haltmeier M 2018 Inverse Problems 34 024006
  * [35] Sandbichler M, Krahmer F, Berer T, Burgholzer P and Haltmeier M 2015 SIAM Journal on Applied Mathematics 75 2475–2494
  * [36] Haltmeier M, Berer T, Moon S and Burgholzer P 2016 Journal of Optics 18 114004
  * [37] Haltmeier M, Neumann L and Rabanser S 2015 Inverse Problems 31 065005
  * [38] Dean-Ben X L, Buehler A, Ntziachristos V and Razansky D 2012 IEEE transactions on medical imaging 31 1922–1928
  * [39] Gutta S, Kalva S K, Pramanik M and Yalavarthy P K 2018 Medical physics 45 3749–3767
  * [40] Prakash J, Sanny D, Kalva S K, Pramanik M and Yalavarthy P K 2018 IEEE transactions on medical imaging 38 1935–1947
  * [41] Bhatt M, Gutta S and Yalavarthy P K 2016 JOSA A 33 1785–1792
  * [42] Provost J and Lesage F 2008 IEEE transactions on medical imaging 28 585–594
  * [43] Betcke M M, Cox B T, Huynh N, Zhang E Z, Beard P C and Arridge S R 2017 IEEE Transactions on Computational Imaging 3 710–721
  * [44] Guo Z, Li C, Song L and Wang L V 2010 Journal of biomedical optics 15 021311
  * [45] Rudin L I, Osher S and Fatemi E 1992 Physica D: nonlinear phenomena 60 259–268
  * [46] Arridge S, Beard P, Betcke M, Cox B, Huynh N, Lucka F, Ogunlade O and Zhang E 2016 Physics in Medicine & Biology 61 8908
  * [47] Han Y, Tzoumas S, Nunes A, Ntziachristos V and Rosenthal A 2015 Medical physics 42 5444–5452
  * [48] Wang K, Su R, Oraevsky A A and Anastasio M A 2012 Physics in Medicine & Biology 57 5399
  * [49] Rosenthal A, Jetzfellner T, Razansky D and Ntziachristos V 2012 IEEE transactions on medical imaging 31 1346–1357
  * [50] Zhang Y, Wang Y and Zhang C 2012 Ultrasonics 52 1046–1055
  * [51] Bredies K, Kunisch K and Pock T 2010 SIAM Journal on Imaging Sciences 3 492–526
  * [52] Hammernik K, Pock T and Nuster R 2017 Variational photoacoustic image reconstruction with spatially resolved projection data Photons Plus Ultrasound: Imaging and Sensing 2017 vol 10064 (International Society for Optics and Photonics) p 100643I
  * [53] Rejesh N, Kalva S, Pramanik M and Arigovindan M 2020 Journal of Instrumentation 15 P12028
  * [54] Prakash J, Raju A S, Shaw C B, Pramanik M and Yalavarthy P K 2014 Biomedical optics express 5 1363–1377
  * [55] Prakash J, Mandal S, Razansky D and Ntziachristos V 2019 IEEE Transactions on Biomedical Engineering 66 2604–2616
  * [56] Blu T and Luisier F 2007 IEEE Transactions on Image Processing 16 2778–2786
  * [57] Xue F, Luisier F and Blu T 2013 IEEE Transactions on Image Processing 22 1954–1968
  * [58] Li J, Luisier F and Blu T 2018 IEEE Transactions on Image Processing 27 92–105
  * [59] Ramani S, Blu T and Unser M 2008 IEEE Transactions on image processing 17 1540–1554
  * [60] Xue F, Yagola A G, Liu J and Meng G 2016 Inverse Problems in Science and Engineering 24 625–646 (Preprint https://doi.org/10.1080/17415977.2015.1054822) URL https://doi.org/10.1080/17415977.2015.1054822
  * [61] Wang K, Ermilov S A, Su R, Brecht H P, Oraevsky A A and Anastasio M A 2011 IEEE transactions on medical imaging 30 203–214
  * [62] Rosenthal A, Razansky D and Ntziachristos V 2010 IEEE transactions on medical imaging 29 1275–1285
  * [63] Tao M, Yang J and He B 2009 TR0918, Department of Mathematics, Nanjing University
  * [64] Wang Y, Yang J, Yin W and Zhang Y 2008 SIAM Journal on Imaging Sciences 1 248–272
  * [65] Arigovindan M, Fung J C, Elnatan D, Mennella V, Chan Y H M, Pollard M, Branlund E, Sedat J W and Agard D A 2013 Proceedings of the National Academy of Sciences 201315675
  * [66] Bertsekas D 1999 Nonlinear Programming (Athena Scientific)
  * [67] Beck A 2017 First-Order Methods in Optimization (MOS-SIAM series on Optimization) chap 6, pp 129–177
  * [68] Liu H, Wang K, Peng D, Li H, Zhu Y, Zhang S, Liu M and Tian J 2016 IEEE transactions on medical imaging 35 2546–2557