# Graph Laplacian for image deblurring111This is a preprint.

Davide Bianchi  Dipartimento di Scienze e Alta Tecnologia
Università dell’Insubria
via Valleggio 11
I-22100 Como, Italy<EMAIL_ADDRESS>,  Alessandro Buccini
Department of Mathematics and Computer Science
University of Cagliari
09124 Cagliari, Italy<EMAIL_ADDRESS>,  Marco Donatelli
Dipartimento di Scienze e Alta Tecnologia
Università dell’Insubria
via Valleggio 11
I-22100 Como, Italy<EMAIL_ADDRESS>and  Emma Randazzo
Dipartimento di Scienze e Alta Tecnologia
Università dell’Insubria
via Valleggio 11
I-22100 Como, Italy<EMAIL_ADDRESS>
###### Abstract.

Image deblurring is relevant in many fields of science and engineering. To
solve this problem, many different approaches have been proposed and among the
various methods, variational ones are extremely popular. These approaches are
characterized by substituting the original problem with a minimization one
where the functional is composed of two terms, a data fidelity term and a
regularization term. In this paper we propose, in the classical
$\ell^{2}-\ell^{1}$ minimization with the non-negativity constraint of the
solution, the use of the graph Laplacian as regularization operator. Firstly,
we describe how to construct the graph Laplacian from the observed noisy and
blurred image. Once the graph Laplacian has been built, we solve efficiently
the proposed minimization problem splitting the convolution operator and the
graph Laplacian by the alternating direction method of multipliers (ADMM).
Some selected numerical examples show the good performances of the proposed
algorithm.

##  1\. Introduction

We are concerned with the solution of the image deblurring problem; see, e.g.,
[31] for more details on image deblurring. We will assume that the blurring is
space-invariant obtaining a linear system of equations

$A\mathbf{x}=\mathbf{b},$

where $\mathbf{x}\in{\mathbb{R}}^{N}$ and $\mathbf{b}\in{\mathbb{R}}^{M}$ are
samplings of the unknown image to recover and the blurred image, respectively,
while $A\in{\mathbb{R}}^{M\times N}$ is a structured matrix (see below) whose
singular values decay rapidly and with no significant gap. The discretization
process, along with measurement errors, introduces some perturbations in the
data, namely $\boldsymbol{\eta}\in{\mathbb{R}}^{M}$ such that
$\|\mathbf{b}-\boldsymbol{\eta}\|=\delta$ in the Euclidean norm, leading to
the system

$A\mathbf{x}=\mathbf{b}+\boldsymbol{\eta}=\mathbf{b}^{\delta}.$

The perturbations $\boldsymbol{\eta}$ are often referred to as noise. Since,
in general, $\boldsymbol{\eta}\notin\mathcal{R}(A)$, we would like to solve
the least-square problem

(1.1) $\arg\min_{\mathbf{x}}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|,$

where $\left\|\cdot\right\|$ is the Euclidean norm. Let $A^{\dagger}$ denote
the Moore-Penrose pseudo-inverse of $A$. The naive solution of (1.1),
$A^{\dagger}\mathbf{b}^{\delta}$, is usually a poor approximation of the
desired solution $\mathbf{x}^{\dagger}=A^{\dagger}\mathbf{b}$; see, e.g., [24,
30] for more details. This is due to the fact that $A$ is severely ill-
conditioned and the solution of (1.1) is extremely sensitive to the presence
of the noise in the data. To compute an approximation of
$\mathbf{x}^{\dagger}$ we need to use regularization methods. These methods
aim at reducing the sensitivity mentioned before. One of the most popular
regularization method is Tikhonov regularization where the original
minimization problem is substituted by another one of the form

(1.2)
$\arg\min_{\mathbf{x}}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|^{2}+\mu\left\|L\mathbf{x}\right\|^{2},$

where $L\in{\mathbb{R}}^{p\times n}$ is an operator such that
$\mathcal{N}(A)\cap\mathcal{N}(L)=\\{\mathbf{0}\\}$. The matrix $L$ is the so-
called regularization operator and its role is to enforce some a-priori
knowledge on the reconstruction. If $A$ is the discretization of an integral
operator, then usually $L$ is chosen to be a discretization of either the
first or the second derivative; see, e.g., [23]. The formulation (1.2) can be
extended for a general $\ell_{p}$-norm

(1.3)
$\arg\min_{\mathbf{x}}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|^{2}+\mu\left\|L\mathbf{x}\right\|^{p}_{p},$

where $\left\|\mathbf{x}\right\|_{p}^{p}=\sum_{i=1}^{n}|x_{i}|^{p}$ for $p>0$.
Note that, for $p<1$ the function
$\mathbf{x}\mapsto\left\|\mathbf{x}\right\|_{p}$ is not a norm; see, e.g.,
[32, 36, 14, 16, 12, 25, 21]. In this paper, we consider le graph Laplacian
for $L$ in (1.3) with $p=1$. Therefore, our minimization problem is of the
form

(1.4) $\arg\min_{\mathbf{x}\geq
0}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|^{2}+\mu\left\|L\mathbf{x}\right\|_{1},$

where we have introduced the non-negativity constraint since images cannot
attain negative values.

Recently, the graph Laplacian of a network, built from the image itself, has
been proposed as regularization operator mainly for image denoising; see.
e.g., [39, 34, 40, 35, 37, 44, 42]. In this paper, we build such an operator
and use it in (1.4) to reconstruct blurred images. We propose an automated
procedure for the construction of an appropriate graph Laplacian and show its
performances in some selected numerical examples. We compare the graph
Laplacian with the standard Total Variation (TV) operator (see [41]) and show
that our proposal can lead to substantial improvements in the quality of the
reconstructed images.

This paper is structured as follows: in Section 2 we recall the definition of
the Laplacian of a given graph and we construct the one that we use in the
following. Section 3 presents our algorithmic proposal for the solution of
(1.4) and Section 4 contains some selected numerical experiments. Finally, we
draw some conclusions in Section 5.

###  1.1. Notation

Discretized images are made by union of several pixels in the plane, and
therefore are well represented by nonnegative two-dimensional discrete
functions $\mathbf{x}:{\mathbb{R}}^{n_{1}\times n_{2}}\to[0,+\infty)$,
$\mathbf{x}(i_{1},i_{2})=x_{i_{1},i_{2}}\in{\mathbb{R}}_{+}$ for
$i_{1}=1,\ldots,n_{1},\,i_{2}=1,\ldots,n_{2}$. Clearly, this choice is not
unique. Since all the operations and analysis we will produce are invariant
with respect to the ordering, for the sake of notational simplicity, we fix
$n_{1}=n_{2}=n=\sqrt{N}$ and consider the lexicographic one-dimensional
ordering, that is, $(i_{1},i_{2})=:i<j:=(j_{1},j_{2})$ if $i_{1}<j_{1}$ or
$i_{1}=j_{1}$ and $i_{2}<j_{2}$. With this choice, an image reads as
$\mathbf{x}:{\mathbb{R}}^{N}\to[0,+\infty)$,
$\mathbf{x}(i)=x_{i}\in{\mathbb{R}}_{+}$ for $i=1,\ldots,n$, and it is said
that $\mathbf{x}$ is the vectorization of a square image.

##  2\. Construction of the graph Laplacian

In this section, we first describe how to construct the Laplacian of a given
weighted graph. Then we show how to build an appropriate graph, i.e., a graph
whose Laplacian is a “good” regularization operator, given a good
approximation of the exact image $\mathbf{x}^{\dagger}$. Finally, we provide
an algorithm to construct our $L$ given the problem (1.1).

Given a countable measure space $\left(\mathcal{V},\nu\right)$, where $\nu$ is
a positive measure, a symmetric non-negative function
$\omega:\mathcal{V}\times\mathcal{V}\to[0,+\infty)$ with zero diagonal is
called undirected _graph_ on $\mathcal{V}$. The elements $i,j$ of the set
$\mathcal{V}$ are called _vertices_ , and two vertices are _connected_ if
$\omega(i,j)>0$. The positive value $w(i,j)$ is called _weight_ associated to
the _edge_ $\\{i,j\\}$; for a modern analytic introduction to graph theory we
refer to [33]. If $\mathcal{V}$ is a finite set of $n$ elements, then the
graph $\omega$ can be uniquely represented, unless permutations, by the
_adjacency matrix_ $\Omega\in{\mathbb{R}}^{N\times N}$

$\left(\Omega\right)_{i,j}:=\omega(i,j).$

The linear operator $L_{\omega}:C(\mathcal{V})\to C(\mathcal{V})$, acting on
the space
$C(\mathcal{V}):=\left\\{\mathbf{x}:\mathcal{V}\to{\mathbb{R}}\right\\}\simeq{\mathbb{R}}^{N}$
via

$L_{\omega}\mathbf{x}(i):=\frac{1}{\nu(i)}\sum_{j}\omega(i,j)\left(\mathbf{x}(i)-\mathbf{x}(j)\right)$

is the graph Laplacian on $C(\mathcal{V})$ associated to the graph $\omega$.
It is a symmetric operator with respect to the inner product

$\langle\mathbf{x},\mathbf{y}\rangle:=\sum_{i}\mathbf{x}(i)\mathbf{y}(i)\nu(i).$

In many applications, a quite standard choice for the measure $\nu$ is the
_degree_ function $\nu=\deg$, defined by

$\deg(i):=\sum_{j}\omega(i,j).$

It measures the whole intensity of the weights associated to the vertex $i$.
Clearly, this choice makes $L_{\omega}$ not symmetric with respect to the
standard Euclidean inner product. A good compromise is to choose the
homogeneous measure associated to the Frobenius norm of $\Omega$, i.e.,
$\nu(i)\equiv\|\Omega\|_{F}$. Let us observe that, writing $D$ as the diagonal
matrix such that $(D)_{i,i}=\deg(i)$, then it is easy to check that

$\|D\|_{1}\leq\|\Omega\|_{F}\leq\|D\|_{2}.$

Henceforth, we will assume $\nu=\|\Omega\|_{F}$. In matrix form, then the
graph Laplacian reads

$L_{\omega}=\frac{D-\Omega}{\|\Omega\|_{F}}.$

We wish to construct a graph $\omega$ so that $L_{\omega}$ can be used in
(1.4). In principle we would like to construct $\omega$ such that

$\left\|L_{\omega}\mathbf{x}^{\dagger}\right\|\approx 0.$

To this aim let $\mathbf{x}^{*}$ be a good approximation of
$\mathbf{x}^{\dagger}$. Define $\omega$ as the weighted and undirected graph
on $\mathcal{V}$, whose nodes are the pixels of $\mathbf{x}^{\dagger}$ and the
weights are defined by

$\omega(i,j)=\left\\{\begin{array}[]{ll}{\rm
e}^{-(x^{*}(i)-x^{*}(j))^{2}/\sigma}&\mbox{if }i\neq j\mbox{ and
}\left\|i-j\right\|_{\infty}\leq R,\\\ 0&\mbox{otherwise,}\end{array}\right.$

where $\sigma>0$ and $R\in{\mathbb{N}}$ are user-defined parameters. Let us
recall that we are using the one-dimensional lexicographic ordering, but the
nodes of $\mathcal{V}$ represent points in ${\mathbb{R}}^{2}$. Therefore,
$i=(i_{1},i_{2}),j=(j_{1},j_{2})$ and
$\|i-j\|_{\infty}=\max\left\\{|i_{1}-j_{1}|;|i_{2}-j_{2}|\right\\}$.
Intuitively, the graph is constructed as follows: we connect two pixels if
they are close enough and we weight their connection depending on how similar
their values are. In particular, we give a strong connection to pixels that
have similar values. The parameter $R$ determines how large is the
neighborhood we consider for each pixel and $\sigma$ determines how strong the
connection between two close pixels should be.

The construction of this graph, and consequently of the graph Laplacian, in
turn depends on the construction of an appropriate approximation
$\mathbf{x}^{*}$. As we show in Section 4, if we could choose
$\mathbf{x}^{*}=\mathbf{x}^{\dagger}$ we would obtain an almost optimal
result. However, this is not possible in realistic scenarios. Therefore, we
wish to provide a practical way to determine a good enough $\mathbf{x}^{*}$ in
a totally automatic way. To compute $\mathbf{x}^{*}$ we propose to solve (1.2)
with a regularization operator defined as follows. Let
$L_{1}\in{\mathbb{R}}^{n\times n}$ be

$L_{1}=\begin{bmatrix}-1&1\\\ &-1&1\\\ &&\ddots&\ddots\\\ &&&-1&1\\\
1&&&&-1\end{bmatrix},$

i.e., $L_{1}$ is a discretization of the first derivative with periodic
boundary conditions (BCs). Let $I_{n}$ be the identity matrix of order $n$, we
define $L_{\rm TV}$ by

(2.1) $L_{\rm TV}=\begin{bmatrix}L_{1}\otimes I_{n}\\\ I_{n}\otimes
L_{1}\end{bmatrix}\in{\mathbb{R}}^{2N\times N}.$

Note that $L_{\rm TV}$ is an extremely sparse matrix formed by two Block
Circulant with Circulant Blocks (BCCB) matrices stacked one over the other.
Therefore, matrix-vector products involving $L_{\rm TV}$ can be performed
extremely cheaply (in particular, the flop count is $O(N)$) and $L^{T}_{\rm
TV}L_{\rm TV}$ is a BCCB matrix. We exploit the latter property below.

To simplify the computations we impose periodic BCs to the matrix $A$. Thanks
to this choice, $A$ is a BCCB matrix; see [31] for more details. We recall
that BCCB matrices are diagonalized by the two-dimensional Fourier matrix. Let
$F_{1}\in{\mathbb{R}}^{n\times n}$ be the Fourier matrix, i.e.,
$(F_{1})_{j,k}={\rm e}^{2\pi\iota(j-1)(k-1)/n}$ with $\iota^{2}=-1$, then the
two-dimensional Fourier matrix is defined by $F=F_{1}\otimes F_{1}$. Note that
matrix-vector products with $F$ and its inverse $F^{*}$ can be performed in
$O(N\log N)$ flops with the aid of the fft and ifft algorithms.

As discussed above we wish to solve (1.2) with $L$ described above to
determine $L_{\omega}$, i.e., we wish to solve a problem of the form

(2.2)
$\mathbf{x}_{\mu}=\arg\min_{\mathbf{x}}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|^{2}+\mu\left\|L_{\rm
TV}\mathbf{x}\right\|^{2},$

for a certain $\mu>0$. Thanks to the structure of $A$ and $L_{\rm TV}$ this
can be solved cheaply for any $\mu$. We can write

(2.3) $A=F^{*}\Sigma F\quad\mbox{and}\quad L_{\rm
TV}=\begin{bmatrix}F^{*}\Lambda_{x}F\\\ F^{*}\Lambda_{y}F\end{bmatrix},$

where $\Sigma$, $\Lambda_{x}$, and $\Lambda_{y}$ are diagonal matrices whose
diagonal entries are the eigenvalues of $A$, $L_{1}\otimes I_{n}$, and
$I\otimes L_{1}$, respectively. We recall that the eigenvalue of a BCCB matrix
$C$ can be computed by $F\mathbf{c}_{1}$, where $\mathbf{c}_{1}$ is the first
column of $C$. Assuming that $\mathcal{N}(A)\cap\mathcal{N}(L_{\rm
TV})=\\{\mathbf{0}\\}$ we have that

$\displaystyle\mathbf{x}_{\mu}$ $\displaystyle=(A^{T}A+\mu L^{T}_{\rm
TV}L_{\rm TV})^{-1}A^{T}\mathbf{b}^{\delta}$
$\displaystyle=\left(F^{*}\Sigma^{*}FF^{*}\Sigma
F+\mu\begin{bmatrix}F^{*}\Lambda_{x}^{*}F&F^{*}\Lambda_{y}^{*}F\end{bmatrix}\begin{bmatrix}F^{*}\Lambda_{x}F\\\
F^{*}\Lambda_{y}F\end{bmatrix}\right)^{-1}F^{*}\Sigma^{*}F\mathbf{b}^{\delta}$
$\displaystyle=\left(F^{*}\Sigma^{*}\Sigma F+\mu
F^{*}(\Lambda_{x}^{*}\Lambda_{x}+\Lambda_{y}^{*}\Lambda_{y})F\right)^{-1}F^{*}\Sigma^{*}F\mathbf{b}^{\delta}$
$\displaystyle=F^{*}\left(\Sigma^{*}\Sigma+\mu(\Lambda_{x}^{*}\Lambda_{x}+\Lambda_{y}^{*}\Lambda_{y})\right)^{-1}\Sigma^{*}F\mathbf{b}^{\delta},$

where the matrix to be inverted is a diagonal matrix. Therefore,
$\mathbf{x}_{\mu}$ can be computed for any $\mu$ cheaply.

We now wish to determine in an automatic way the parameter $\mu$. We employ
the Generalized Cross Validation (GCV).

Denote by $G(\mu)$ the following function

$G(\mu)=\frac{\left\|A\mathbf{x}_{\mu}-\mathbf{b}^{\delta}\right\|^{2}}{\mbox{{trace}}(I-A(A^{T}A+\mu
L^{T}_{\rm TV}L_{\rm TV})^{-1}A^{T})^{2}}.$

The GCV parameter is the minimizer of $G(\mu)$, i.e., $\mu_{\rm
GCV}=\arg\min_{\mu}G(\mu)$. Given the decomposition (2.3) the value of
$G(\mu)$ can be computed in a straightforward way. Introduce the following
notation

$r_{\mu}=\left\|A\mathbf{x}_{\mu}-\mathbf{b}^{\delta}\right\|\quad\mbox{and}\quad
t_{\mu}=\mbox{{trace}}(I-A(A^{T}A+\mu L^{T}_{\rm TV}L_{\rm TV})^{-1}A^{T}),$

therefore, $G(\mu)=r_{\mu}^{2}/t_{\mu}^{2}$. Using the spectral decomposition
of $A$ we have

$\displaystyle r_{\mu}$
$\displaystyle=\left\|A\mathbf{x}_{\mu}-\mathbf{b}^{\delta}\right\|=\left\|F^{*}\Sigma
FF^{*}\left(\Sigma^{*}\Sigma+\mu(\Lambda_{x}^{*}\Lambda_{x}+\Lambda_{y}^{*}\Lambda_{y})\right)^{-1}\Sigma^{*}F\mathbf{b}^{\delta}-\mathbf{b}^{\delta}\right\|$
$\displaystyle=\left\|\Sigma\left(\Sigma^{*}\Sigma+\mu(\Lambda_{x}^{*}\Lambda_{x}+\Lambda_{y}^{*}\Lambda_{y})\right)^{-1}\Sigma^{*}F\mathbf{b}^{\delta}-F\mathbf{b}^{\delta}\right\|$
$\displaystyle=\left\|(\Sigma\Sigma^{*}\left(\Sigma^{*}\Sigma+\mu(\Lambda_{x}^{*}\Lambda_{x}+\Lambda_{y}^{*}\Lambda_{y})\right)^{-1}-I)\widehat{\mathbf{b}^{\delta}}\right\|,$

where $\widehat{\mathbf{b}^{\delta}}=F\mathbf{b}^{\delta}$. We now move to the
computation of $t_{\mu}$

$\displaystyle t_{\mu}$ $\displaystyle=\mbox{{trace}}(I-A(A^{T}A+\mu
L^{T}_{\rm TV}L_{\rm TV})^{-1}A^{T})$
$\displaystyle=\mbox{{trace}}(I-F^{*}\Sigma(\Sigma^{*}\Sigma+\mu\left(\Lambda_{x}^{*}\Lambda_{x}+\Lambda_{y}^{*}\Lambda_{y})\right)^{-1}\Sigma^{*}F)$
$\displaystyle=\mbox{{trace}}(I-\Sigma\Sigma^{*}(\Sigma^{*}\Sigma+\mu\left(\Lambda_{x}^{*}\Lambda_{x}+\Lambda_{y}^{*}\Lambda_{y})\right)^{-1}).$

We can observe that, once the decompositions (2.3) and
$\widehat{\mathbf{b}^{\delta}}$ have been computed, the evaluation of $G(\mu)$
can be done in $O(N)$ flops. This allows for an extremely fast determination
of $\mu_{\rm GCV}$. Finally, we select as $\mathbf{x}^{*}$ the solution of the
minimization problem

(2.4)
$\mathbf{x}^{*}=\arg\min_{\mathbf{x}}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|^{2}+\mu_{\rm
GCV}\left\|L_{\rm TV}\mathbf{x}\right\|^{2}.$

######  Remark 2.1.

If $A$ is constructed with BCs different from the periodic ones it is still
possible to compute a fairly accurate approximation of $G(\mu)$ using Krylov
subspace methods; see [26, 27, 17, 8].

We summarize the procedure to construct $L_{\omega}$ in Algorithm 1.

Algorithm 1 Construction of $L_{\omega}$ for image deblurring

  * 1:

Input: $A\in{\mathbb{R}}^{N\times N}$,
$\mathbf{b}^{\delta}\in{\mathbb{R}}^{N}$, $R>0$, $\sigma>0$

  * 2:

Output: $L_{\omega}\in{\mathbb{R}}^{n\times n}$

  * 3:

Construct $L_{\rm TV}$ as defined in (2.1)

  * 4:

$\Sigma=\mbox{{diag}}\left(F(A_{(:,1)})\right)$

  * 5:

$\Lambda_{x}=\mbox{{ diag}}\left(F((L_{\rm TV})_{(:,1:N)})\right)$

  * 6:

$\Lambda_{y}=\mbox{{ diag}}\left(F((L_{\rm TV})_{(:,N+1:2N)})\right)$

  * 7:

$\widehat{\mathbf{b}^{\delta}}=F\mathbf{b}^{\delta}$

  * 8:

$\mu_{\rm GCV}=\arg\min_{\mu}\frac{r_{\mu}^{2}}{t_{\mu}^{2}}$, where

$\left\\{\begin{array}[]{l}r_{\mu}=\left\|(\Sigma\Sigma^{*}\left(\Sigma^{*}\Sigma+\mu(\Lambda_{x}^{*}\Lambda_{x}+\Lambda_{y}^{*}\Lambda_{y})\right)^{-1}-I)\widehat{\mathbf{b}^{\delta}}\right\|\\\
t_{\mu}=\mbox{{trace}}(I-\Sigma\Sigma^{*}(\Sigma^{*}\Sigma+\mu\left(\Lambda_{x}^{*}\Lambda_{x}+\Lambda_{y}^{*}\Lambda_{y})\right)^{-1})\end{array}\right.$

  * 9:

$\mathbf{x}^{*}=F^{*}\left(\Sigma^{*}\Sigma+\mu(\Lambda_{x}^{*}\Lambda_{x}+\Lambda_{y}^{*}\Lambda_{y})\right)^{-1}\Sigma^{*}\widehat{\mathbf{b}^{\delta}}$

  * 10:

Construct $\Omega\in{\mathbb{R}}^{n\times n}$ as

$\omega(i,j)=\left\\{\begin{array}[]{ll}{\rm
e}^{-(x^{*}(i)-x^{*}(j))^{2}/\sigma}&\mbox{if }i\neq j\mbox{ and
}\left\|i-j\right\|_{\infty}\leq R,\\\ 0&\mbox{else,}\end{array}\right.$

with $i$ and $j$ representing two-dimensional indexes, sorted in lexicographic
order

  * 11:

$D=\mbox{{diag}}\\{\sum_{j=1}^{n}(\Omega)_{(i,j)}\\}$

  * 12:

$L_{\omega}=\frac{D-\Omega}{\left\|\Omega\right\|_{F}}$

##  3\. Graph Laplacian deblurring

We now describe the non-linear model we employ to compute an approximate
solution of (1.1). We consider the graph Laplacian $L_{\omega}$ constructed by
Algorithm 1 and we use it in (1.4). Therefore, we wish to solve the following

(3.1) $\arg\min_{\mathbf{x}\geq
0}\frac{1}{2}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|^{2}+\mu\left\|L_{\omega}\mathbf{x}\right\|_{1}.$

To solve this problem we use the Alternating Direction Multiplier Method
(ADMM); see, e.g., [6] for a recent review. We use ADMM since it allows us to
decouple the $\ell^{2}$ and $\ell^{1}$ norms as well as the matrices $A$ and
$L_{\omega}$. The latter point is extremely relevant since, as we discuss
below, both matrices have exploitable structures, however, they are difficult
to exploit together.

We first reformulate (3.1) in an equivalent form

$\arg\min_{\mathbf{x},\mathbf{y},\mathbf{w},\mathbf{z}}\left\\{\frac{1}{2}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|^{2}+\mu\left\|\mathbf{z}\right\|_{1}+\iota_{0}(\mathbf{w}),\mbox{
s.t.
}\mathbf{x}=\mathbf{y},\mathbf{x}=\mathbf{w},\mathbf{z}=L_{\omega}\mathbf{y},\right\\}$

or equivalently

(3.2)
$\begin{split}\arg\min_{\mathbf{x},\mathbf{y},\mathbf{w},\mathbf{z}}&\frac{1}{2}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|^{2}+\mu\left\|\mathbf{z}\right\|_{1}+\iota_{0}(\mathbf{w})\\\
\mbox{ s.t. }&\begin{bmatrix}I&O\\\ O&I\\\
I&O\end{bmatrix}\begin{bmatrix}\mathbf{x}\\\
\mathbf{z}\end{bmatrix}-\begin{bmatrix}I&O\\\ L_{\omega}&O\\\
O&I\end{bmatrix}\begin{bmatrix}\mathbf{y}\\\
\mathbf{w}\end{bmatrix}=\mathbf{0},\end{split}$

where $O$ and $\mathbf{0}$ denote the zero matrix and the zero vector,
respectively, and $\iota_{0}$ is the indicator function of the nonnegative
cone, i.e.,

$\iota_{0}(\mathbf{x})=\left\\{\begin{array}[]{ll}0&\mbox{if }\mathbf{x}\geq
0,\\\ +\infty&\mbox{otherwise.}\end{array}\right.$

We can construct the augmented Lagrangian of (3.2)

$\displaystyle\mathcal{L}_{\rho}(\mathbf{x},\mathbf{y},\mathbf{w},\mathbf{z};\boldsymbol{\lambda})=$
$\displaystyle\frac{1}{2}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|^{2}+\mu\left\|\mathbf{z}\right\|_{1}+\iota_{0}(\mathbf{w})$
$\displaystyle-\left\langle\boldsymbol{\lambda},\begin{bmatrix}I&O\\\ O&I\\\
I&O\end{bmatrix}\begin{bmatrix}\mathbf{x}\\\
\mathbf{z}\end{bmatrix}-\begin{bmatrix}I&O\\\ L_{\omega}&O\\\
O&I\end{bmatrix}\begin{bmatrix}\mathbf{y}\\\
\mathbf{w}\end{bmatrix}\right\rangle$
$\displaystyle+\frac{\rho}{2}\left\|\begin{bmatrix}I&O\\\ O&I\\\
I&O\end{bmatrix}\begin{bmatrix}\mathbf{x}\\\
\mathbf{z}\end{bmatrix}-\begin{bmatrix}I&O\\\ L_{\omega}&O\\\
O&I\end{bmatrix}\begin{bmatrix}\mathbf{y}\\\
\mathbf{w}\end{bmatrix}\right\|^{2},$

where $\rho>0$ is a fixed parameter and
$\boldsymbol{\lambda}\in{\mathbb{R}}^{3N}$ is the Lagrangian multiplier.
Applying ADMM we get the iterations

$\left\\{\begin{array}[]{l}\vspace*{2mm}\begin{bmatrix}\mathbf{x}^{(k+1)}\\\
\mathbf{z}^{(k+1)}\end{bmatrix}=\arg\min_{\mathbf{x},\mathbf{z}}\mathcal{L}_{\rho}(\mathbf{x},\mathbf{y}^{(k)},\mathbf{w}^{(k)},\mathbf{z};\boldsymbol{\lambda}^{(k)}),\\\
\begin{bmatrix}\mathbf{y}^{(k+1)}\\\
\mathbf{w}^{(k+1)}\end{bmatrix}=\arg\min_{\mathbf{y},\mathbf{w}}\mathcal{L}_{\rho}(\mathbf{x}^{(k+1)},\mathbf{y},\mathbf{w},\mathbf{z}^{(k+1)};\boldsymbol{\lambda}^{(k)}),\\\
\boldsymbol{\lambda}^{(k+1)}=\boldsymbol{\lambda}^{(k)}+\rho\left(\begin{bmatrix}I&O\\\
O&I\\\ I&O\end{bmatrix}\begin{bmatrix}\mathbf{x}^{(k+1)}\\\
\mathbf{z}^{(k+1)}\end{bmatrix}-\begin{bmatrix}I&O\\\ L_{\omega}&O\\\
O&I\end{bmatrix}\begin{bmatrix}\mathbf{y}^{(k+1)}\\\
\mathbf{w}^{(k+1)}\end{bmatrix}\right).\end{array}\right.$

We can write $\boldsymbol{\lambda}=\begin{bmatrix}\boldsymbol{\lambda}_{1}\\\
\boldsymbol{\lambda}_{2}\\\ \boldsymbol{\lambda}_{3}\end{bmatrix}$ with
$\boldsymbol{\lambda}_{j}\in{\mathbb{R}}^{N}$ for $j=1,2,3$. Therefore the
minimization problems above decouples and we obtain

$\left\\{\begin{array}[]{l}\mathbf{x}^{(k+1)}=\arg\min_{\mathbf{x}}\frac{1}{2}\left\|A\mathbf{x}-\mathbf{b}^{\delta}\right\|^{2}-\left\langle\begin{bmatrix}\boldsymbol{\lambda}_{1}^{(k)}\\\
\boldsymbol{\lambda}^{(k)}_{3}\end{bmatrix},\begin{bmatrix}\mathbf{x}-\mathbf{y}^{(k)}\\\
\mathbf{x}-\mathbf{w}^{(k)}\end{bmatrix}\right\rangle+\frac{\rho}{2}\left\|\begin{bmatrix}\mathbf{x}-\mathbf{y}^{(k)}\\\
\mathbf{x}-\mathbf{w}^{(k)}\end{bmatrix}\right\|^{2},\\\
\mathbf{z}^{(k+1)}=\arg\min_{\mathbf{z}}\mu\left\|\mathbf{z}\right\|_{1}-\left\langle\boldsymbol{\lambda}_{2}^{(k)},\mathbf{z}-L_{\omega}\mathbf{y}^{(k)}\right\rangle+\rho_{2}\left\|\mathbf{z}-L_{\omega}\mathbf{y}^{(k)}\right\|,\\\
\mathbf{y}^{(k+1)}=\arg\min_{\mathbf{y}}\frac{\rho}{2}\left\|\mathbf{x}^{(k+1)}-\mathbf{y}\right\|^{2}-\left\langle\begin{bmatrix}\boldsymbol{\lambda}_{1}^{(k)}\\\
\boldsymbol{\lambda}_{2}^{(k)}\end{bmatrix},\begin{bmatrix}\mathbf{x}^{(k+1)}-\mathbf{y}\\\
L_{\omega}\mathbf{z}^{(k+1)}-\mathbf{y}\end{bmatrix}\right\rangle+\frac{\rho}{2}\left\|\begin{bmatrix}\mathbf{x}^{(k+1)}-\mathbf{y}\\\
L_{\omega}\mathbf{z}^{(k+1)}-\mathbf{y}\end{bmatrix}\right\|^{2},\\\
\mathbf{w}^{(k+1)}=\arg\min_{\mathbf{w}}\frac{\rho}{2}\left\|\mathbf{x}^{(k+1)}-\mathbf{w}\right\|^{2}-\left\langle\boldsymbol{\lambda}_{3}^{(k)},\mathbf{x}^{(k+1)}-\mathbf{w}\right\rangle+\iota_{0}(\mathbf{w}).\end{array}\right.$

All of these minimization problems have a closed form for their solution,
namely

$\left\\{\begin{array}[]{l}\mathbf{x}^{(k+1)}=(A^{T}A+2\rho
I)^{-1}(A^{T}\mathbf{b}^{\delta}+\rho\mathbf{y}^{(k)}-\boldsymbol{\lambda}^{(k)}_{1}+\rho\mathbf{w}^{(k)}-\boldsymbol{\lambda}^{(k)}_{3}),\\\
\mathbf{z}^{(k+1)}=S_{\mu/\rho}\left(L_{\omega}\mathbf{y}^{(k)}-\boldsymbol{\lambda}^{(k)}_{2}/\rho\right),\\\
\mathbf{y}^{(k+1)}=(L_{\omega}^{T}L_{\omega}+I)^{-1}(L_{\omega}^{T}(\mathbf{z}^{(k+1)}+\boldsymbol{\lambda}^{(k)}_{2}/\rho)+\mathbf{x}^{(k+1)}+\boldsymbol{\lambda}^{(k)}_{1}/\rho),\\\
\mathbf{w}^{(k+1)}=\left(\mathbf{x}^{(k+1)}+\boldsymbol{\lambda}_{3}/\rho\right)_{+},\\\
\boldsymbol{\lambda}_{1}^{(k+1)}=\boldsymbol{\lambda}_{1}^{(k)}+\rho(\mathbf{x}^{(k+1)}-\mathbf{y}^{(k+1)}),\\\
\boldsymbol{\lambda}_{2}^{(k+1)}=\boldsymbol{\lambda}_{2}^{(k)}+\rho(\mathbf{z}^{(k+1)}-L_{\omega}\mathbf{y}^{(k+1)}),\\\
\boldsymbol{\lambda}_{3}^{(k+1)}=\boldsymbol{\lambda}_{3}^{(k)}+\rho(\mathbf{x}^{(k+1)}-\mathbf{w}^{(k+1)}),\end{array}\right.$

where $S_{\mu}$ denotes the soft-thresholding operator with parameter $\mu$,
i.e., $S_{\mu}(\mathbf{x})=\mbox{{sign}}(\mathbf{x})(|\mathbf{x}|-\mu)_{+}$,
where the operations are meant element-wise and $(x)_{+}=\max\\{x,0\\}$ is the
metric projection into the nonnegative cone. Note that each iteration requires
the solution of two linear systems. The linear system involving $A$ can be
easily solved using the fft algorithm if periodic BCs are employed; see above.
If other BCs are employed, the structure of the matrix, in general, does not
allow to use fast transform for the solution of the system. Nevertheless, this
system can be solved by an iterative method using a circulant preconditioner.
On the other hand, the solution of linear system with the $L_{\omega}$ matrix
can be easily computed using the lsqr algorithm applied to the equivalent
least-square problem since the matrix $L_{\omega}$ is extremely sparse.

We would like to briefly discuss the use of the lsqr method for the solution
of the system

(3.3)
$(L_{\omega}^{T}L_{\omega}+I)\mathbf{y}^{(k+1)}=(L_{\omega}^{T}(\mathbf{z}^{(k+1)}+\boldsymbol{\lambda}^{(k)}_{2}/\rho)+\mathbf{x}^{(k+1)}+\boldsymbol{\lambda}^{(k)}_{1}/\rho.$

The linear system of equations (3.3) is equivalent to the least squares
problem

(3.4)
$\begin{split}\mathbf{y}^{(k+1)}&=\arg\min_{\mathbf{y}}\left\|\begin{bmatrix}L_{\omega}\\\
I\end{bmatrix}\mathbf{y}-\begin{bmatrix}\mathbf{z}^{(k+1)}+\boldsymbol{\lambda}^{(k)}_{2}/\rho\\\
\mathbf{x}^{(k+1)}+\boldsymbol{\lambda}^{(k)}_{1}/\rho\end{bmatrix}\right\|^{2}\\\
&=\arg\min_{\mathbf{y}}\left\|\widehat{L}\mathbf{y}-\widehat{\mathbf{v}}^{(k)}\right\|^{2}.\end{split}$

The lsqr algorithm is an iterative method that determines an approximate
solution of (3.4) into a Krylov subspace. In particular, at its $j-$th
iteration the lsqr algorithm determines a solution of (3.4) in the Krylov
subspace
$\mathcal{K}_{j}(\widehat{L}^{T}\widehat{L},\widehat{L}^{T}\widehat{\mathbf{v}}^{(k)})$,
where

$\mathcal{K}_{j}(\widehat{L}^{T}\widehat{L},\widehat{L}^{T}\widehat{\mathbf{v}}^{(k)})=\mbox{{span}}\left\\{\widehat{L}^{T}\widehat{\mathbf{v}}^{(k)},(\widehat{L}^{T}\widehat{L})\widehat{L}^{T}\widehat{\mathbf{v}}^{(k)},\ldots,(\widehat{L}^{T}\widehat{L})^{j-1}\widehat{L}^{T}\widehat{\mathbf{v}}^{(k)}\right\\}.$

The lsqr method requires one matrix-vector product with $\widehat{L}$ and one
with $\widehat{L}^{T}$ at each iteration. Therefore, since $\widehat{L}$ is
extremely sparse the flop count per iteration is $O(N)$. Moreover lsqr is
mathematically equivalent to the cg method applied to (3.3). However, its
implementation is more stable. Nevertheless, the number of iteration required
to converge is proportional to $\kappa(\widehat{L})$ which is extremely small;
see, e.g., [5] for a discussion on lsqr and cg.

We summarize our approach in Algorithm 2.

Algorithm 2 Graph Laplacian image deblurring

  * 1:

Input: $A\in{\mathbb{R}}^{N\times N}$,
$\mathbf{b}^{\delta}\in{\mathbb{R}}^{N}$, $R>0$, $\sigma>0$, $\rho>0$,
$\tau>0$, $K>0$, $\mu>0$

  * 2:

Output: $\mathbf{x}\in{\mathbb{R}}^{n\times n}$

  * 3:

Run Algorithm 1 to compute $L_{\omega}$

  * 4:

$\mathbf{y}^{(0)}=\mathbf{w}^{(0)}=\boldsymbol{\lambda}_{1}^{(0)}=\boldsymbol{\lambda}_{2}^{(0)}=\boldsymbol{\lambda}_{3}^{(0)}=\mathbf{0}$

  * 5:

for $k=0,\ldots,K$

  * 6:

$\mathbf{x}^{(k+1)}=(A^{T}A+2\rho
I)^{-1}(A^{T}\mathbf{b}^{\delta}+\rho\mathbf{y}^{(k)}-\boldsymbol{\lambda}^{(k)}_{1}+\rho\mathbf{w}^{(k)}-\boldsymbol{\lambda}^{(k)}_{3})$

  * 7:

$\mathbf{z}^{(k+1)}=S_{\mu/\rho}\left(L_{\omega}\mathbf{y}^{(k)}-\boldsymbol{\lambda}^{(k)}_{2}/\rho\right)$

  * 8:

$\mathbf{y}^{(k+1)}=(L_{\omega}^{T}L_{\omega}+I)^{-1}(L_{\omega}^{T}(\mathbf{z}^{(k+1)}+\boldsymbol{\lambda}^{(k)}_{2}/\rho)+\mathbf{x}^{(k+1)}+\boldsymbol{\lambda}^{(k)}_{1}/\rho)$

  * 9:

$\mathbf{w}^{(k+1)}=P_{0}\left(\mathbf{x}^{(k+1)}+\boldsymbol{\lambda}_{3}/\rho\right)$

  * 10:

$\boldsymbol{\lambda}_{1}^{(k+1)}=\boldsymbol{\lambda}_{1}^{(k)}+\rho(\mathbf{x}^{(k+1)}-\mathbf{y}^{(k+1)})$

  * 11:

$\boldsymbol{\lambda}_{2}^{(k+1)}=\boldsymbol{\lambda}_{2}^{(k)}+\rho(\mathbf{z}^{(k+1)}-L_{\omega}\mathbf{y}^{(k+1)})$

  * 12:

$\boldsymbol{\lambda}_{3}^{(k+1)}=\boldsymbol{\lambda}_{3}^{(k)}+\rho(\mathbf{x}^{(k+1)}-\mathbf{w}^{(k+1)})$

  * 13:

if
$k>1\;\&\;\left\|\mathbf{x}^{k+1}-\mathbf{x}^{k}\right\|\leq\tau\left\|\mathbf{x}^{k}\right\|$

  * 14:

Exit

  * 15:

end if

  * 16:

end for

  * 17:

$\mathbf{x}=\mathbf{x}^{(k+1)}$

Note that since the functional (1.4) is convex we can apply the following
classical result on ADMM

######  Theorem 3.1 (see, e.g., Section 3.2 of Boyd et al. [6]).

With the notation of Algorithm 2 it holds that

  1. (i)

$\lim_{k\rightarrow\infty}\left\|\begin{bmatrix}I&O\\\ O&I\\\
I&O\end{bmatrix}\begin{bmatrix}\mathbf{x}^{(k)}\\\
\mathbf{z}^{(k)}\end{bmatrix}-\begin{bmatrix}I&O\\\ L_{\omega}&O\\\
O&I\end{bmatrix}\begin{bmatrix}\mathbf{y}^{(k)}\\\
\mathbf{w}^{(k)}\end{bmatrix}\right\|=0$, i.e., the iterates approach
feasibility as $k\rightarrow\infty$;

  2. (ii)

$\lim_{k\rightarrow\infty}\frac{1}{2}\left\|A\mathbf{x}^{(k)}-\mathbf{b}^{\delta}\right\|^{2}+\mu\left\|\mathbf{z}^{(k)}\right\|_{1}+\iota_{0}(\mathbf{w}^{(k)})=p^{*}$,
where $p^{*}$ is the minimum of (3.1);

  3. (iii)

$\lim_{k\rightarrow\infty}\boldsymbol{\lambda}_{k}=\boldsymbol{\lambda}^{*}$,
where $\boldsymbol{\lambda}^{*}$ is a dual optimal point, i.e., a saddle point
of $L_{0}$.

######  Remark 3.2.

ADMM can be slow to converge in certain scenarios. It is not in the scope of
this paper to propose a fast algorithm for the solution of (3.1). Rather we
wish to show the potentiality of $L_{\omega}$ as a regularization operator.
Nevertheless, it is possible to accelerate the convergence of ADMM by
extrapolation methods to improve the convergence rate of ADMM; see, e.g., [9,
29].

##  4\. Numerical Examples

We now report some selected numerical examples to show the performances of the
proposed method. We are particularly interested in showing that the graph
Laplacian constructed in Algorithm 2 provides better reconstructions than the
classical TV approach.

We compare the results obtained using $L=L_{\rm TV}$ and $L=L_{\omega}$ in
(1.4) with the solution $\mathbf{x}^{*}$ computed in line $9$ of Algorithm 2.
To compute the solution of (1.4) with $L=L_{\rm TV}$ we use the algorithm
described in [22]. Moreover, we want to show the fill potentiality of the
proposed approach. To this aim we construct the operator $L_{\omega}$ using
the exact image $\mathbf{x}_{\rm true}$ and we denote it by
$\widetilde{L}_{\omega}$. Obviously this approach is not feasible in realistic
scenarios, nevertheless it allows us to show the full potentiality of the
proposed method.

In all our experiments we set the parameters as specified in Table 1. We
compare the considered method in terms of accuracy using the Relative
Restoration Error (RRE) computed as

${\rm RRE}(\mathbf{x})=\frac{\left\|\mathbf{x}-\mathbf{x}_{\rm
true}\right\|_{2}}{\left\|\mathbf{x}_{\rm true}\right\|_{2}},$

the Pick Signal to Noise Ration (PSNR), defined by

${\rm PSNR}(\mathbf{x})=20\log_{10}\left(\frac{Nm}{\left\|\mathbf{x}_{\rm
true}\right\|_{2}}\right),$

where $m$ denotes the maximum value achievable by $\mathbf{x}_{\rm true}$.
Moreover, we consider the Structure SIMilarity index (SSIM), constructed in
[43]. The definition of the SSIM is extremely involved, here we simply recall
that this statistical index measures how structurally similar two images are,
in particular, the higher the SSIM the more similar the images are, and its
highest achievable value is $1$.

Table 1. Setting of the parameters in Algorithm 2. Paramter | Value | Description
---|---|---
$R$ | $10$ | Support of the weight function in the Graph
$\sigma$ | $10^{-2}$ | Variance of the weight function in the Graph
$\rho$ | $10^{-1}$ | Augmentation parameter in ADMM
$\tau$ | $10^{-4}$ | Stopping criterion for ADMM
$K$ | $3000$ | Maximum number of iterations
$\mu$ | Hand-tuned | Regularization parameter in (1.4)

#### Example 1.

Our first example is the atmosphericBlur50 test case of the RestoreTools
toolbox [2]. We report the exact image, the PSF, and the blurred an noisy
image in Fig. 2. The norm of the noise, denoted by $\delta$, that corrupts the
data is approximately $1\%$ of the norm of the exact right-hand side
$\mathbf{b}^{\delta}$.

We report the obtained results with the considered methods in Table 2. We can
observe that $\ell^{2}-\ell^{1}$ methods provide much more accurate results
than the classical Tikhonov method, especially in terms of SSIM. The
reconstruction obtained with $L=\widetilde{L}_{\omega}$, i.e., using the graph
related to the exact image, is extremely accurate. To validate our model we
show in Fig. 1 a visualization of
$\left|\widetilde{L}_{\omega}\mathbf{x}^{\dagger}\right|$. We can observe that
it is extremely sparse. Therefore, we expect our model to provide accurate
reconstructions. However, this approach is not feasible in real scenarios.
Nevertheless, we can also observe that using $L=L_{\omega}$ improves the
quality of the restoration with respect to the classic TV. This is confirmed
by the visual inspection of the reconstructions in Fig. 3. Comparing the
reconstructions obtained by the $\ell^{2}-\ell^{1}$ methods we can observe
that the choice $L=L_{\rm TV}$ leads to more noisy reconstructions than the
ones obtained with the graph Laplacian.

Figure 1. Example 1: Visualization of
$\left|\widetilde{L}_{\omega}\mathbf{x}^{\dagger}\right|$ in the jet colormap.
The color blue represents the $0$s in the image.

(a)

(b)

(c)

Figure 2. Example 1: (a) True image ($256\times 256$ pixels), (b) PSF
($256\times 256$ pixels), (c) Blurred and noisy image ($256\times 256$ pixels
with $\delta\approx 0.01\left\|\mathbf{b}\right\|_{2}$).

(a)

(b)

(c)

Figure 3. Example 1 reconstructions: (a) $\ell^{2}-\ell^{1}$ with $L=L_{\rm
TV}$, (b) $\ell^{2}-\ell^{1}$ with $L=L_{\omega}$, (c) $\ell^{2}-\ell^{1}$
with $L=\widetilde{L}_{\omega}$.

#### Example 2.

For our second example we consider the Hubble image in Fig. 4(a) and we blur
it with the PSF in Fig. 4(b). We then add white Gaussian noise such that
$\delta=0.1\left\|\mathbf{b}\right\|_{2}$ obtaining the blurred and noisy
image in Fig. 4(c).

We compute approximate solutions with the considered algorithms and report the
obtained RRE, PSNR, and SSIM in Table 2. We can observe that our proposal
provides the best reconstruction both in terms of RRE (and, therefore of PSNR)
and SSIM. This is confirmed by the visual inspection of the reconstructions in
Fig. 5. We would like to stress that, similarly as the previous example, the
unconstrained Tikhonov method computes extremely noisy reconstructions.

(a)

(b)

(c)

Figure 4. Example 2: (a) True image ($256\times 256$ pixels), (b) PSF
($9\times 9$ pixels), (c) Blurred and noisy image ($256\times 256$ pixels with
$\delta=0.1\left\|\mathbf{b}\right\|_{2}$).

(a)

(b)

(c)

Figure 5. Example 2 reconstructions: (a) $\ell^{2}-\ell^{1}$ with $L=L_{\rm
TV}$, (b) $\ell^{2}-\ell^{1}$ with $L=L_{\omega}$, (c) $\ell^{2}-\ell^{1}$
with $L=\widetilde{L}_{\omega}$.

#### Example 3.

For our third example we consider the Saturn image in Fig. 6(a). We blur it
with a non-symmetric PSF (see Fig. 6(b)) and add $5\%$ of white Gaussian
noise, i.e., $\delta=0.05\left\|\mathbf{b}\right\|_{2}$ obtaining the image in
Fig. 6(c).

We report in Table 2 the obtained results with the considered algorithms. We
can observe that our proposal provides a very accurate reconstruction in terms
of RRE and PSNR. However, the SSIM of the computed solution is slightly lower
than the one obtained with the standard TV regularization. In Fig. 7 we report
all the computed solution. We would like to observe that the reconstruction
obtained with the classical TV regularization presents a very heavy stair-case
effect while the approximate solution obtained with our proposal avoids this
issue. To show this we propose in Fig. 8 blow-ups of the central part of the
image of the exact solution and of the reconstructions obtained by TV
regularization and our approach. We can observe that the TV reconstruction
presents a very heavy stair-case effect that is avoided by our proposal. The
solution computed by our method manages to remain extremely sharp, while
avoiding any stair-casing.

(a)

(b)

(c)

Figure 6. Example 3: (a) True image ($256\times 256$ pixels), (b) PSF
($17\times 17$ pixels), (c) Blurred and noisy image ($256\times 256$ pixels
with $\delta=0.05\left\|\mathbf{b}\right\|_{2}$).

(a)

(b)

(c)

Figure 7. Example 3 reconstructions: (a) $\ell^{2}-\ell^{1}$ with $L=L_{\rm
TV}$, (b) $\ell^{2}-\ell^{1}$ with $L=L_{\omega}$, (c) $\ell^{2}-\ell^{1}$
with $L=\widetilde{L}_{\omega}$.

(a)

(b)

(c)

Figure 8. Example 3 blow ups of the exact solution and of two reconstructions
in the jet colormap: (a) True solution, (b) $\ell^{2}-\ell^{1}$ with $L=L_{\rm
TV}$, (c) $\ell^{2}-\ell^{1}$ with $L=L_{\omega}$.

#### Example 4.

We consider the exact image in Fig. 9(a) and we blur it with an average PSF
(see Fig. 9(b)), we add white Gaussian noise so that
$\left\|\mathbf{b}-\mathbf{b}^{\delta}\right\|=0.03\left\|\mathbf{b}\right\|$
obtaining Fig. 9(c). We construct the graph Laplacian related to the exact
solution $\widetilde{L}_{\omega}$ and deblur the blurred and noisy image with
Tikhonov (with $\mu=\mu_{\rm GCV}$) and by minimizing (1.4) with both
$L=L_{\rm TV}$ and $L=\widetilde{L}_{\omega}$. We report the results obtained
in Table 2. We can observe that the proposed algorithm with the perfect choice
of $\widetilde{L}_{\omega}$ largely outperforms the other approaches
furnishing a very accurate reconstruction of the proposed image. Moreover, the
choice of $L=L_{\omega}$, which we recall that it does not require any a
priori information on the exact solution, is still more accurate than the
classical TV. This is confirmed by the visual inspection of the
reconstructions in Fig. 10.

(a)

(b)

(c)

Figure 9. Example 4: (a) True image ($256\times 256$ pixels), (b) PSF
($12\times 12$ pixels), (c) Blurred and noisy image ($256\times 256$ pixels
with $\delta=0.03\left\|\mathbf{b}\right\|_{2}$).

(a)

(b)

(c)

Figure 10. Example 4 reconstructions: (a) $\ell^{2}-\ell^{1}$ with $L=L_{\rm TV}$, (b) $\ell^{2}-\ell^{1}$ with $L=L_{\omega}$, (c) $\ell^{2}-\ell^{1}$ with $L=\widetilde{L}_{\omega}$. Table 2. Comparison of the RRE, PSNR, and SSIM. Example | Method | RRE | PSNR | SSIM
---|---|---|---|---
Example 1 | Tikhonov | $0.22299$ | $26.663$ | $0.55512$
$\ell^{2}-\ell^{1}$ with $L=L_{\rm TV}$ | $0.19152$ | $27.984$ | $0.92623$
$\ell^{2}-\ell^{1}$ with $L=L_{\omega}$ | ${0.17763}$ | ${28.638}$ | ${0.93971}$
$\ell^{2}-\ell^{1}$ with $L=\widetilde{L}_{\omega}$ | $0.083333$ | $35.212$ | $0.98129$
Example 2 | Tikhonov | $0.17352$ | $25.735$ | $0.55241$
$\ell^{2}-\ell^{1}$ with $L=L_{\rm TV}$ | $0.15492$ | $26.720$ | $0.80458$
$\ell^{2}-\ell^{1}$ with $L=L_{\omega}$ | $0.14968$ | $27.019$ | $0.81256$
$\ell^{2}-\ell^{1}$ with $L=\widetilde{L}_{\omega}$ | $0.10096$ | $30.439$ | $0.89943$
Example 3 | Tikhonov | $0.080283$ | $33.715$ | $0.72254$
$\ell^{2}-\ell^{1}$ with $L=L_{\rm TV}$ | $0.068917$ | $35.041$ | $0.94873$
$\ell^{2}-\ell^{1}$ with $L=L_{\omega}$ | $0.060094$ | $36.231$ | $0.94809$
$\ell^{2}-\ell^{1}$ with $L=\widetilde{L}_{\omega}$ | $0.040382$ | $39.684$ | $0.96489$
Example 4 | Tikhonov | $0.16236$ | $27.160$ | $0.73224$
$\ell^{2}-\ell^{1}$ with $L=L_{\rm TV}$ | $0.15299$ | $27.686$ | $0.86899$
$\ell^{2}-\ell^{1}$ with $L=L_{\omega}$ | $0.14716$ | $28.024$ | $0.85765$
$\ell^{2}-\ell^{1}$ with $L=\widetilde{L}_{\omega}$ | $0.085936$ | $32.695$ | $0.93887$

##  5\. Conclusions

In this paper we have proposed a new regularization operator for
$\ell^{2}-\ell^{1}$ minimization. The construction of this operator is
automatic and extremely cheap to perform. We have shown that the proposed
method outperforms the classical TV approach. Matter of future research
include the application of the proposed method to more general inverse
problems as well as the integration of the considered method with the
$\ell^{p}-\ell^{q}$ minimization proposed in [32, 36, 14, 15, 16, 12] or with
iterative regularization methods like, e.g., Linearized Bregman splitting [18,
19, 11, 20, 13] and Iterated Tikhonov with general penalty terms [3, 10, 7,
4]. Another line of future research is the construction of more sophisticated
graphs $\omega$ which can better exploit the structure of the given image
itself. Such constructions may stem from a PDEs approach; see, e.g., [28, 38,
1].

## Acknowledgments

A.B., and M.D. are members of the GNCS-INdAM group. D.B. is member of the
GNAMPA-INdAM group. A.B. research is partially supported by the Regione
Autonoma della Sardegna research project “Algorithms and Models for Imaging
Science [AMIS]” (RASSR57257, intervento finanziato con risorse FSC 2014-2020 -
Patto per lo Sviluppo della Regione Sardegna).

## References

  * [1] Adriani, A., Bianchi, D., Serra-Capizzano, S.: Asymptotic Spectra of Large (Grid) Graphs with a Uniform Local Structure (Part I): Theory.  Milan J. Math. 87(2), 169–199 (2020)
  * [2] Berisha, S., Nagy, J.G.: Iterative methods for image restoration.  Tech. rep., Department of Mahtematics and Computer Science, Emory University (2012).  http://www.mathcs.emory.edu/~nagy/RestoreTools/IR.pdf
  * [3] Bianchi, D., Buccini, A., Donatelli, M., Serra-Capizzano, S.: Iterated fractional Tikhonov regularization.  Inverse Problems 31(5), 055005 (2015)
  * [4] Bianchi, D., Donatelli, M.: On generalized iterated Tikhonov regularization with operator-dependent seminorms.  Electron. Trans. Numer. Anal. 47, 73–99 (2017)
  * [5] Björck, Å.: Numerical methods for least squares problems.  SIAM (1996)
  * [6] Boyd, S., Parikh, N., Chu, E.: Distributed optimization and statistical learning via the alternating direction method of multipliers.  Now Publishers Inc (2011)
  * [7] Buccini, A.: Regularizing preconditioners by non-stationary iterated tikhonov with general penalty term.  Applied Numerical Mathematics 116, 64–81 (2017)
  * [8] Buccini, A.: Generalized Cross Validation stopping rule for Iterated Tikhonov regularization.  In preparation (2021)
  * [9] Buccini, A., Dell’Acqua, P., Donatelli, M.: A general framework for admm acceleration.  Numerical Algorithms 85, 829–848 (2020)
  * [10] Buccini, A., Donatelli, M., Reichel, L.: Iterated Tikhonov with general penalty term.  Accepted on: Numerical Linear Algebra and Applications (2016)
  * [11] Buccini, A., Park, Y., Reichel, L.: Numerical aspects of the nonstationary modified linearized Bregman algorithm.  Applied Mathematics and Computation 337, 386–398 (2018)
  * [12] Buccini, A., Pasha, M., Reichel, L.: Modulus-based iterative methods for constrained $\ell^{p}-\ell^{q}$ minimization.  Inverse Problems 36(8), 084001 (2020)
  * [13] Buccini, A., Pasha, M., Reichel, L.: Linearized Krylov subspace Bregman iteration with nonnegativity constraint.  Numerical Algorithms (In press)
  * [14] Buccini, A., Reichel, L.: An $\ell^{2}-\ell^{q}$ regularization method for large discrete ill-posed problems.  Journal of Scientific Computing 78(3), 1526–1549 (2019)
  * [15] Buccini, A., Reichel, L.: An $\ell^{2}$-$\ell^{q}$ regularization method for large discrete ill-posed problems.  Journal of Scientific Computing 78(3), 1526–1549 (2019)
  * [16] Buccini, A., Reichel, L.: An $\ell^{p}$-$\ell^{q}$ minimization method with cross-validation for the restoration of impulse noise contaminated images.  Journal of Computational and Applied Mathematics 375, 112824 (2020)
  * [17] Buccini, A., Reichel, L.: Generalized cross validation for $\ell^{p}$-$\ell^{q}$ minimization.  Under Review (2021)
  * [18] Cai, J.F., Osher, S., Shen, Z.: Linearized Bregman iterations for frame-based image deblurring.  SIAM Journal on Imaging Sciences 2(1), 226–252 (2009)
  * [19] Cai, J.F., Osher, S., Shen, Z.: Split Bregman methods and frame based image restoration.  Multiscale Modeling & Simulation 8(2), 337–369 (2009)
  * [20] Cai, Y., Donatelli, M., Bianchi, D., Huang, T.Z.: Regularization preconditioners for frame-based image deblurring with reduced boundary artifacts.  SIAM Journal on Scientific Computing 38(1), B164–B189 (2016)
  * [21] Chan, R.H., Liang, H.X.: Half-quadratic algorithm for $\ell_{p}$-$\ell_{q}$ problems with applications to TV-$\ell_{1}$ image restoration and compressive sensing.  In: Efficient Algorithms for Global Optimization Methods in Computer Vision, pp. 78–103. Springer (2014)
  * [22] Chan, R.H., Tao, M., Yuan, X.: Constrained total variation deblurring models and fast algorithms based on alternating direction method of multipliers.  SIAM Journal on imaging Sciences 6(1), 680–697 (2013)
  * [23] Donatelli, M., Reichel, L.: Square smoothing regularization matrices with accurate boundary conditions.  Journal of Computational and Applied Mathematics 272, 334–349 (2014)
  * [24] Engl, H.W., Hanke, M., Neubauer, A.: Regularization of Inverse Problems.  Kluwer, Doordrecht (1996)
  * [25] Estatico, C., Gratton, S., Lenti, F., Titley-Peloquin, D.: A conjugate gradient like method for p-norm minimization in functional spaces.  Numerische Mathematik 137(4), 895–922 (2017)
  * [26] Fenu, C., Reichel, L., Rodriguez, G.: GCV for Tikhonov regularization via global Golub-Kahan decomposition.  Numerical Linear Algebra with Applications 23(3), 467–484 (2016)
  * [27] Fenu, C., Reichel, L., Rodriguez, G., Sadok, H.: GCV for Tikhonov regularization by partial SVD.  BIT Numerical Mathematics 57, 1019–1039 (2017)
  * [28] Gilboa, G., Osher, S.: On generalized iterated Tikhonov regularization with operator-dependent seminorms.  SIAM Multiscale Model. Simul. 7(3), 1005–1028 (2009)
  * [29] Goldstein, T., O’Donoghue, B., Setzer, S., Baraniuk, R.: Fast alternating direction optimization methods.  SIAM Journal on Imaging Sciences 7(3), 1588–1623 (2014)
  * [30] Hansen, P.C.: Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion.  SIAM (1998)
  * [31] Hansen, P.C., Nagy, J.G., O’Leary, D.P.: Deblurring Images: Matrices, Spectra, and Filtering.  SIAM, Philadelphia (2006)
  * [32] Huang, G., Lanza, A., Morigi, S., Reichel, L., Sgallari, F.: Majorization-minimization generalized Krylov subspace methods for $\ell_{p}$-$\ell_{q}$ optimization applied to image restoration.  BIT Numerical Mathematics pp. 1–28 (2017)
  * [33] Keller, M., Lenz, D., Wojciechowski, R.K.: Graphs and discrete Dirichlet spaces.  Grundlehren der mathematischen Wissenschaften: Springer, forthcoming. (2021)
  * [34] Kheradmand, A., Milanfar, P.: Motion deblurring with graph laplacian regularization.  In: Digital Photography XI, vol. 9404, p. 94040C. International Society for Optics and Photonics (2015)
  * [35] Krishnan, D., Fergus, R.: Fast image deconvolution using hyper-laplacian priors.  In: Advances in neural information processing systems, pp. 1033–1041 (2009)
  * [36] Lanza, A., Morigi, S., Reichel, L., Sgallari, F.: A generalized Krylov subspace method for $\ell_{p}$-$\ell_{q}$ minimization.  SIAM Journal on Scientific Computing 37(5), S30–S50 (2015)
  * [37] Li, F., Ng, M.K.: Image colorization by using graph bi-laplacian.  Advances in Computational Mathematics 45(3), 1521–1549 (2019)
  * [38] Lou, Y., Zhang, X., Osher, S., Bertozzi, A.: Image recovery via nonlocal operators.  J. Sci. Comput. 42, 185–197 (2010)
  * [39] Meyer, F.G., Shen, X.: Perturbation of the eigenvectors of the graph laplacian: Application to image denoising.  Applied and Computational Harmonic Analysis 36(2), 326–334 (2014)
  * [40] Pang, J., Cheung, G.: Graph laplacian regularization for image denoising: Analysis in the continuous domain.  IEEE Transactions on Image Processing 26(4), 1770–1785 (2017)
  * [41] Rudin, L.I., Osher, S., Fatemi, E.: Nonlinear total variation based noise removal algorithms.  Physica D: nonlinear phenomena 60(1-4), 259–268 (1992)
  * [42] Susnjara, A., Perraudin, N., Kressner, D., Vandergheynst, P.: Accelerated filtering on graphs using lanczos method.  arXiv preprint arXiv:1509.04537 (2015)
  * [43] Wang, Z., Bovik, A.C., Sheikh, H.R., Simoncelli, E.P.: Image quality assessment: from error visibility to structural similarity.  IEEE Transactions on Image Processing 13(4), 600–612 (2004)
  * [44] Yağan, A.C., Özgen, M.T.: A spectral graph wiener filter in graph fourier domain for improved image denoising.  In: 2016 IEEE Global Conference on Signal and Information Processing (GlobalSIP), pp. 450–454. IEEE (2016)


# Absolute Properties of the Oscillating Eclipsing Algol X Trianguli

Jae Woo Lee1, Kyeongsoo Hong2, Jang-Ho Park1, Marek Wolf3, and Dong-Jin Kim1
1Korea Astronomy and Space Science Institute, Daejeon 34055, Republic of Korea
2Institute for Astrophysics, Chungbuk National University, Cheongju 28644,
Republic of Korea  3Astronomical Institute, Faculty of Mathematics and
Physics, Charles University in Prague, 180 00 Praha 8, V Holešovičkách 2,
Czech Republic<EMAIL_ADDRESS>
###### Abstract

We report results from the TESS photometric data and new high-resolution
spectra of the Algol system X Tri showing short-period pulsations. From the
echelle spectra, the radial velocities of the eclipsing pair were measured,
and the rotational rate and effective temperature of the primary star were
obtained to be $v_{1}$$\sin$$i=84\pm 6$ km s-1 and $T_{\rm eff,1}=7900\pm 110$
K, respectively. The synthetic modeling of these observations implies that X
Tri is in synchronous rotation and is physically linked to a visual companion
TIC 28391715 at a separation of about 6.5 arcsec. The absolute parameters of
our target star were accurately and directly determined to be $M_{1}=2.137\pm
0.018$ $M_{\odot}$, $M_{2}=1.101\pm 0.010$ $M_{\odot}$, $R_{1}=1.664\pm 0.010$
$R_{\odot}$, $R_{2}=1.972\pm 0.010$ $R_{\odot}$, $L_{1}=9.67\pm 0.55$
$L_{\odot}$, and $L_{2}=2.16\pm 0.09$ $L_{\odot}$. The phase-binned mean light
curve was used to remove the binary effect from the observed TESS data.
Multifrequency analysis of the residuals revealed 16 significant frequencies,
of which the high-frequency signals between 37 day-1 and 48 day-1 can be
considered probable pulsation modes. Their oscillation periods of
0.021$-$0.027 days and pulsation constants of 0.014$-$0.018 days are typical
values of $\delta$ Sct variables. The overall results demonstrate that X Tri
is an oEA star system, consisting of a $\delta$ Sct primary and its lobe-
filling companion in the semi-detached configuration.

##  1 INTRODUCTION

Most stars in binary and multiple systems are born almost simultaneously and
evolved by interaction between their components (Duchêne & Kraus 2013). The
pulsations in such systems are very useful for probing the impacts of mass
transfer and tides on stellar evolution and interior structure (Murphy 2018;
Bowman et al. 2019). Many of them are $\delta$ Sct pulsations in the mass-
accreting primary components of semi-detached interacting Algols (Zhou 2010;
Lee et al. 2016; Kahraman Aliçavuş et al. 2017; Liakos & Niarchos 2017;
Mkrtichian et al. 2018), the so-called oscillating eclipsing Algol (oEA) class
(Mkrtichian et al. 2002, 2004). In oEA precursors, the original more massive
components fill their Roche lobes and some masses are transferred to the less
massive companions. As a consequence of the mass transfer and accretion, the
gainers become the present A/F-type pulsating primary components and the
losers become the current low-mass secondary stars of spectral type F$-$K that
are magnetically active.

For eclipsing $\delta$ Sct stars, a possible relationship between the orbital
and pulsation periods was established empirically and theoretically,
respectively, by Soydugan et al. (2006) and Zhang et al. (2013); the longer
the former ($P_{\rm orb}$), the longer the latter ($P_{\rm pul}$). Liakos &
Niarchos (2015, 2017) suggested that the oscillations are influenced by the
binarity in $P_{\rm orb}<$ 13 days. However, Kahraman Aliçavuş et al. (2017)
showed that the threshold is too low and the $\delta$ Sct pulsators in
eclipsing binaries (EBs) oscillate in shorter pulsation periods with lower
amplitudes than single variables. Mkrtichian et al. (2018) reported that the
rapid mass accretion on the oscillating primary in the oEA system RZ Cas is
produced by the magnetic cycle of its lobe-filling companion, and this can
induce a periodic variation in pulsation properties. Recently, many new
$\delta$ Sct-like pulsators in EBs have been detected from the archival data
of the TESS mission, but their physical properties remain unclear (e.g.,
Kahraman Aliçavuş et al. 2022; Shi et al. 2022). A review of such pulsating
EBs was presented in Lampens (2021).

We have been carrying out follow-up spectroscopy of the EB systems with
pulsating components. The Algol system X Tri (TIC 28391714; Gaia EDR3
298268990328628224; TYC 1763-2733-1; $T_{\rm p}$ = $+$8.565; $V\rm_{T}$ =
$+$9.036, $(B-V)\rm_{T}$ = $+$0.328; $P_{\rm orb}$ = 0.9715 days) has been
included as a probable oscillating EB in our target list. It was individually
discovered by Walker (1921) and Neujmin (1922) as a variable star. The
previous spectra of the program target were secured by Struve (1946). He
classified the spectral types of both components as A3$+$G3 and obtained only
primary radial velocities (RVs) with a semi-amplitude of $K_{1}$ = 110 km s-1
and a rotation rate of $v_{1}$ = 50 km s-1, about 0.56 times slower than its
synchronous value (van Hamme & Wilson 1990). Bozkurt et al. (1976) and Liakos
et al. (2010) made complete photoelectric and CCD light curves in the $BV$ and
$BVRI$ bands, respectively, and analyzed their own observations. Mezzetti et
al. (1980) solved again the multiband light curves of Bozkurt et al. (1976)
with the elements of spectroscopic orbits of Struve (1946) and presented the
masses and radii of each component: $M_{1}$ = 2.3 $\pm$ 0.7 $M_{\odot}$,
$M_{2}$ = 1.2 $\pm$ 0.3 $M_{\odot}$, $R_{1}$ = 1.71 $\pm$ 0.03 $R_{\odot}$,
and $R_{2}$ = 1.96 $\pm$ 0.03 $R_{\odot}$. The historical results indicate
that the program target has a semi-detached configuration with the lobe-
filling secondary. Meanwhile, its orbital period has been studied by Gadomski
(1932), Wood (1950), Wood & Forbes (1963), Frieboes-Conde & Herczeg (1973),
Rafert (1982), Rovithis-Livaniou et al. (2000), and Qian (2002). Most
recently, Liakos et al. (2010) suggested that X Tri is likely a quintuple
system with three circumbinary companions in periods and minimal masses of
36.9 yrs and 0.18 M⊙, 22.4 yrs and 0.24 M⊙, and 16.8 yrs and 0.22 M⊙.

There have been some photometric attempts to search for pulsation features in
X Tri. Kim et al. (2003) and Liakos & Niarchos (2009) did not find any
noticeable signals from their 2001 and 2007$-$2008 observations, respectively,
while Turner et al. (2014) detected an oscillation with a period of 0.0220
days and an amplitude of 20 mmag in $V$ passband during their campaign on
October 2010. This article is the tenth in a series studying the absolute
properties of the pulsating EBs by combining our high-resolution spectra with
existing or new photometric data (cf. Hong et al. 2015; Kim et al. 2022).

##  2 TESS PHOTOMETRY AND ECLIPSE TIMES

X Tri was observed in Sector 17 (BJD 2,458,764.688 $-$ 2,458,789.684) by
camera 1 of TESS (Ricker et al. 2015) and recorded with a 2-min cadence. The
observations had an interruption of about 1.5 days near the midpoint for data
download. For this study, we used the SAP-FLUX data available from the MAST
portal111https://archive.stsci.edu/ and did not include all points between BJD
2,458,772.466 and 2,458,776.278 showing unrealistic light variations. The raw
data were detrended and normalized by applying a second-order polynomial fit
to each of the two segments divided by the interruption (Lee et al. 2019). We
transformed the corrected fluxes to magnitudes.

From the time-series data, 38 mid-eclipse timings were calculated with the
Kwee & van Woerden (1956) method and are compiled in Table 1. With the timing
measurements, we derived the linear eclipse ephemeris for the TESS
observations, as follows:

$\mbox{Min I}=\mbox{BJD}~{}2,458,765.148638(5)+0.9715257(3)E.$ (1)

The orbital period ($P_{\rm orb}$) is equivalent to a frequency of $f_{\rm
orb}$ = 1.0293088 $\pm$ 0.0000005 day-1. Using this ephemeris, the phase-
folded light curve of X Tri is presented in Figure 1.

Using all available minimum epochs, we plotted the current eclipse timing
$O-C$ diagram of X Tri, which is given in the upper panel of Figure 2. In this
study we collected 653 times of minimum light. Most of the eclipse times can
be found f.e. in the $O-C$ Gateway
database222http://var2.astro.cz/ocgate/index.php?lang=en. The general trend of
the $O-C$ values is a downward parabola, suggesting a long-term period
decrease in this system. With the least-squares method, we derived the
following quadratic ephemeris:

$C=\mbox{HJD}~{}2,447,296.2617(3)+0.97152997(5)E-1.91(2)\times 10^{-10}E^{2}.$
(2)

As was claimed by previous investigators, the secular period decrease cannot
be connected with a mass transfer in the system. In case of a conservative
mass transfer in such a semidetached system with the Roche-lobe filling
secondary and more massive primary component, the orbital period should be
increasing. The timing residuals from the quadratic ephemeris are presented in
the lower panel. One can see irregular or quasi-periodic oscillations with
typical lengths of about 20 years, which cannot be explained simply by a
light-time effect (LITE) due to a circumbinary object. In this context, Liakos
et al. (2010) suggested a multiple LITE of three additional bodies, which
partially explain these variations.

X Tri is a member of a similar short-period Algol group (TU Her, FH Ori, TY
Peg, Z Per, Y Psc) with a known long-term period decrease as well as short-
period changes (Qian 2002). The nature of the long-term decrease as well as
short-term irregular period changes has not yet been clearly and convincingly
explained, and are often attributed to the mass and angular momentum loss
mechanism, as developed by Tout & Hall (1991). Concerning a possible magnetic
activity cycle in the secondary companion, there is not enough energy in the
outer layers of the active cool star to reproduce the observed orbital period
changes (Rovithis-Livaniou et al. 2000).

##  3 NEW SPECTROSCOPY AND DATA ANALYSIS

We conducted time-series spectroscopy through the BOES spectrograph mounted on
the 1.8 m telescope at the Bohyunsan Optical Astronomy Observatory (Kim et al.
2007). Thirty-five spectra of X Tri were acquired on five nights between
December 2019 and November 2020. Their wavelength range and resolving power
are $\lambda$ = 3600$-$10,200 $\rm\AA$ and $R$ = 30,000, respectively. The
integration time of the program target was set to 30 min, corresponding to
0.021 of the eclipse period. The raw CCD images were processed utilizing the
IRAF packages CCDRED and ECHELLE. The reduction process was identical to that
employed by Hong et al. (2015). The consequential signal-to-noise ratio (SNR)
is approximately 40, measured between 4000 $\rm\AA$ and 5000 $\rm\AA$.

Struve (1946) reported that his results were not fully exempted from
rotational broadening and line blending between the components. Trailed
spectra are of great use to reveal and track the orbital motions of binary
stars (e.g., Lee et al. 2018, 2020a,b). In order to look for sets of
absorption lines for X Tri, we made the phase-folded trailed spectra from the
BOES observations and examined them in detail at full wavelength. As a
consequence, we found that Fe I $\lambda$4957.61 and Ca I doublet
$\lambda$6103/6122 lines are strong enough to individually trace both
component stars.

In Figure 3, the trailed spectra of the Ca I doublet are displayed in the
upper panel and the spectrum observed on BJD 2,459,169.0452 (0.73 phase) is
presented as an example in the lower panel. For RV determinations, the
absorption lines were fitted several times by means of double Gaussian
functions with the deblending routine in the IRAF $splot$ task. The final RVs
and their 1-$\sigma$ values were averaged and standardized from these
measurements. They are presented in Table 2 and Figure 4, where triangles and
circles denote Fe I and Ca I RVs, respectively. From a sine-curve fit to each
star’s RVs, we obtained the preliminary semi-amplitudes of $K_{1}$ = 107.8 km
s-1 and $K_{2}$ = 202.2 km s-1, respectively, corresponding to a mass ratio of
$q$ = 0.533.

The atmospheric parameters of X Tri can be measured by comparing its observed
spectra and synthetic models (Hong et al. 2017; Lee et al. 2022b). First of
all, we constructed a disentangled spectrum from our BOES spectra and RV
measures using the FDBinary code of Ilijić et al. (2004). Then, the projected
rotational velocity ($v_{1}$$\sin$$i$) and effective temperature ($T_{\rm
eff,1}$) of the primary component were yielded by minimizing the $\chi^{2}$
statistic between the reconstructed spectrum and model spectra. To do this, we
generated grids of synthetic spectra in ranges of 10 $\leq$ $v_{1}$$\sin$$i$
$\leq$ 200 km s-1 and 6000 $\leq$ $T_{\rm eff,1}$ $\leq$ 10,000 K from the
BOSZ spectral library (Bohlin et al. 2017). The surface gravity of $\log
g_{1}$ = 4.3 (cf. Section 6) and the solar metallicity of $[$Fe/H$]=0$ were
applied to the calculations of the model spectra, because the two parameters
are highly correlated with stellar temperatures. For the $\chi^{2}$
minimization, we selected seven spectral regions (Ca II $\lambda$3933, Hδ
$\lambda$4101, Hγ $\lambda$4340, Fe I $\lambda$4383, Hβ $\lambda$4861, Ca I
doublet $\lambda$6103/6122, Hα $\lambda$6563), which are displayed in Figure
5. We acquired the best-fitting parameters of $v_{1}$$\sin$$i=84\pm 6$ km s-1
and $T_{\rm eff,1}=7900\pm 110$ K. These values indicate that the primary star
is rotating $\sim$1.7 times faster than Struve’s (1946) rate of 50 km s-1 and
its spectral type is deduced to be A6.5 V (Pecaut & Mamajek 2013).

On the other hand, there is no evidence of any Balmer emission nor the
presence of a circumstellar matter in our observed spectra. This may be
because the emission feature is very weak and/or our SNR is not enough to
detect it.

##  4 BINARY MODELING AND ABSOLUTE DIMENSIONS

Like the historical data (Bozkurt et al. 1976; Liakos et al. 2010), the TESS
light curve of X Tri has similarities to that of classical Algols. To see if
there was any light variability in X Tri, we divided the time-series TESS data
into 19 light curves at intervals of $P_{\rm orb}$ = 0.9715257 days and
measured each of them at four orbital phases 0.0 (Min I), 0.25 (Max I), 0.5
(Min II), and 0.75 (Max II). These measurements appear in Table 3, where the
BJD time of each segment is the middle of the beginning and the end. The
average values at each phase are given in the last line of Table 3. As
presented in this table, the primary eclipse depth is more than five times
($\sim$1.24 mag) deeper than the secondary one and the depth difference
implies a large temperature difference between the component stars. Further,
Max I is $\sim$0.006 mag brighter than Max II, which may be indicative of
stellar activities on and around both components. The light levels for each
light curve were subtracted from the TESS average for the whole datasets and
the mean residuals are plotted in Figure 6. We show that the primary minima
gradually brightened up after the interruption in the middle point.

The TESS archive data and our velocity curves of X Tri were solved applying
the Wilson-Devinney (W-D) code (Wilson & Devinney 1971, van Hamme & Wilson
2007) with the proximity and eclipse effects on RVs, as in the case of V404
Lyr (Lee et al. 2020a) and WASP 0131+28 (Lee et al. 2020b). The main
difference between the 2007 version we used and the latest version is (1)
direct distance estimation and (2) the use of absolute physical units (e.g.,
Kallrath, J. 2022), which does not affect our X Tri parameters. Thus, we used
the 2007-version W-D code with features optimized by us, including repetitive
routines. The absolute dimensions and distance determination of X Tri were
calculated independently regardless of the W-D code.

In our synthesis, we fixed the primary’s temperature at $T_{\rm eff,1}=7900\pm
110$ K measured from the BOES spectra. It is assumed that X Tri A and B have
predominately radiative and convective atmospheres, respectively. Bolometric
albedo and gravity-darkening exponent were set to be ($A_{1}$, $g_{1}$) =
(1.0, 1.0) and ($A_{2}$, $g_{2}$) = (0.5, 0.32). The logarithmic law was used
for the limb-darkening effect (van Hamme 1993). Based on $v_{1}$$\sin$$i=84\pm
6$ km s-1 and for the lobe-filling secondary, we adopted the rotation
parameters of $F_{1,2}$ = 1.0 in the W-D synthetic code.

In a simultaneous light and RV solution, it is hard to assign reasonable
weights to these observations. To resolve this, the binary modeling in this
article was conducted in two steps. To begin with, the TESS light curve was
analyzed using the $q$ value taken from our RV fitting. In the second, our RV
curves were fitted by adopting the light curve parameters derived in the first
step. Then again, the TESS data were analyzed applying the newly calculated
spectroscopic elements ($q$, $\gamma$, and $a$). This process was carried out
for two cases without and with a starspot, and was repeated until both curves
were reasonably solved.

Table 4 presents the final model parameters with a cool starspot located on
the secondary’s surface, and the RV synthetic curves from them are shown as
red lines in Figure 4. The red line in the top panel of Figure 1 presents our
spot solution, and the light curve residuals corresponding to the unspotted
and cool-spot models are given in the second and third panels. Although the
light diminution seen between phases 0.5 and 0.8 was better described by the
spot model, we could not fit the quasi-sinusoidal modulations to occur around
both eclipses. The ‘W’-shape variation pattern of the residuals could result
from some combination of stellar phenomena, such as magnetic activity, a gas
stream from the lobe-filling component, and its impact on the detached
companion. However, at present, we do not have a clear explanation for this
puzzle. As in our previous studies (Lee et al. 2021, 2022a), the uncertainties
of the binary star parameters in Table 4 were computed following the procedure
applied by Southworth et al. (2020) for the TESS time-series data.

The binary star model indicates that X Tri is a semi-detached Algol, in which
the detached primary component fills about 69 % of its inner Roche lobe. The
absolute dimensions of each star in Table 5 were computed from the synthetic
solution of the TESS data and our double-lined RVs. We used the solar
magnitude of $M_{\rm bol}$⊙ = +4.73 and the empirical bolometric corrections
(BCs) presented by Torres (2010). For measuring the distance to X Tri, we
first obtained its apparent magnitude and color index (CI) to be $V$ = +9.01
$\pm$ 0.02 and ($B-V$) = 0.28 $\pm$ 0.03 mag, respectively, from the Tycho
magnitudes of $B\rm_{T}$ and $V\rm_{T}$ (Høg et al. 2000). Then, the intrinsic
CI of ($B-V$)0 = +0.18 $\pm$ 0.02 was estimated applying $T_{\rm eff,1}$ to
the temperature$-$color calibration (Flower 1996), and $E$($B-V$) = +0.10
$\pm$ 0.04. Our computed values lead to a nominal distance of 215 $\pm$ 6 pc,
which is in good agreement with 210 $\pm$ 2 pc inverted from GAIA EDR3 ($\pi$
= 4.753 $\pm$ 0.042 mas; Gaia Collaboration et al. 2021).

##  5 PULSATIONAL CHARACTERISTICS

The fundamental parameters of X Tri in Table 5 indicate that the hotter
primary star is present in the main-sequence (MS) $\delta$ Sct region of the
Hertzsprung-Russell diagram, while the less massive secondary lies above the
MS band where the lobe-filling companions of other semi-detached Algols exist
(İbanoǧlu et al. 2006; Lee et al. 2016). Turner et al. (2014) found a
pulsation period of 0.0220 $\pm$0.0002 days in our target during their
photometric survey to find eclipsing $\delta$ Sct stars. As displayed in the
third panel of Figure 1, even our spot model did not fully explain the TESS
light curve of X Tri with exceptional precision. It is not easy to extract
reliable pulsating signals from the ‘W’-shape residuals. For useful
multifrequency analyses, we made 1000 normal points phase-binning at intervals
of 0.001. The mean light curve and the residuals from it are presented as the
cyan dots in the top and bottom panels of Figure 1. The corresponding
residuals distributed in BJD are given in Figure 7. We recognized that the
binary signals were satisfactorily removed from the observed data and the
primary eclipse depth varied with time.

In order to investigate the oscillating features of X Tri, the PERIOD04
software by Lenz & Breger (2005) was introduced into the out-of-eclipse light
curve residuals. The amplitude spectrum of the program target was calculated
to the Nyquist frequency limit of $f_{\rm Nyq}\simeq$ 360 day-1. Using the
pre-whitening technique (Lee et al. 2014), we found 16 significant signals
with SNR larger than about 4 for each peak (Breger et al. 1993). The resultant
frequency parameters are given in Table 6 and their uncertainties were
obtained following Kallinger et al. (2008). The synthetic curve prepared from
them appears as a solid line in Figure 7, and the periodogram of X Tri are
presented in Figure 8. There are no conspicuous signals between 60 day-1 and
$f_{\rm Nyq}$.

As illustrated in Table 6 and Figure 8, the main signals of X Tri are
distributed in the two frequency domains of $<$ 3 day-1 and 36$-$52 day-1. Of
our 16 frequencies, the highest amplitude signal $f_{2}$ corresponds to a
frequency of 45.45 $\pm$ 0.41 day-1 detected by Turner et al. (2014) within
their errors, and $f_{14}$ appears to be the orbital harmonic of $2f_{\rm
orb}$. We carefully explored combination frequencies within the Rayleigh
frequency of $R_{\rm f}$ = 0.042 day-1, from which our remarks are given in
the last column of Table 6. We think the low frequencies below $f_{\rm orb}$
might be aliases or artefacts due to systematic trends. In contrast, most high
frequencies may be identified as possible pulsation signals arising from the
hotter and more massive primary.

##  6 DISCUSSION AND CONCLUSIONS

We have analyzed new spectroscopic observations of the semi-detached Algol X
Tri with archival TESS data. The RV curves of the eclipsing pair were obtained
by using the double Gaussian as a fitting function, and we measured the first
RVs for the cool secondary star. The minimum $\chi^{2}$ method was applied to
yield a rotational rate of $v_{1}$$\sin$$i=84\pm 6$ km s-1 and a surface
temperature of $T_{\rm eff,1}=7900\pm 110$ K for the primary star. With the
atmospheric parameters, the TESS light and our RV curves let us characterize
each binary component. From their simultaneous modeling, we determined the
accurate masses and radii with less than 1% uncertainties, as follows:
$M_{1}=2.137\pm 0.018$ $M_{\odot}$, $M_{2}=1.101\pm 0.010$ $M_{\odot}$,
$R_{1}=1.664\pm 0.010$ $R_{\odot}$, and $R_{2}=1.972\pm 0.010$ $R_{\odot}$.
For X Tri A, the synchronous rotation of $v_{\rm 1,sync}$ = 86.7 $\pm$ 0.5 km
s-1 is well matched with our measurement of $v_{1}$$\sin$$i$, unlike Struve’s
(1946). At present, our program target is expected to be in a state of
synchronous rotation. Meanwhile, X Tri is a semidetached Algol in which the
less massive secondary fills the inner critical Roche lobe, while the more
massive primary is detached and its fill-out factor is about 0.69. Thus, our
target X Tri is a regular Algol, not a reverse Algol, i.e. the massive primary
fills its Roche lobe.

Our binary synthesis revealed a third light of $l_{3}$ = 0.0374, which could
originate from either circumbinary objects, as suggested by eclipse timing
analyses, or other stars near X Tri because each pixel size of the TESS CCD
camera is 21 arcsec on the sky (Ricker et al. 2015). There is a neighboring
star TIC 28391715 (2MASS J02003411+2753233; Gaia EDR3 298268986033892480;
$T_{\rm p}$ = $+$12.007) with a separation of about 6.5 arcsec from our target
star and a magnitude difference of 3.442 mag between them in the TESS
passband. The brightness difference indicates that TIC 28391715 is about 24
times fainter than X Tri and hence contributes within 4 % of the total light.
Thus, most of the $l_{3}$ source may be the nearby star TIC 28391715. The
distance to TIC 28391715 is 207.6 $\pm$ 0.7 pc from GAIA EDR3 parallax of
4.816 $\pm$ 0.016 mas (Gaia Collaboration et al. 2021) and this is a good
match to the Gaia distance 210 $\pm$ 2 pc of X Tri. Further, the TESS targets
seem to share a common proper motion; $\mu_{\alpha}\cos\delta$ = 25.893 $\pm$
0.043 mas yr-1 and $\mu_{\delta}$ = $-$14.053 $\pm$ 0.046 mas yr-1 for X Tri,
and $\mu_{\alpha}\cos\delta$ = 26.072 $\pm$ 0.018 mas yr-1 and $\mu_{\delta}$
= $-$12.953 $\pm$ 0.018 mas yr-1 for TIC 28391715. So then, TIC 28391715 may
be gravitationally bound to the X Tri system, and the separation between them
is estimated to be $\sim$1360 AU.

To extract more reliable oscillating frequencies, we formed a normal light
curve at 0.001 phase intervals and the outside-eclipse residuals from it were
applied to the PERIOD04 software package. Sixteen significant frequencies with
SNR $\gtrsim$ 4 were detected from the multifrequency analysis, of which the
low frequencies in the gravity ($g$)-mode domain are thought to be aliases or
artefacts. Except for possible combinations within $R_{\rm f}$, our target
star X Tri pulsates in the pressure ($p$)-mode frequencies between 37 day-1
and 48 day-1 that originate from the primary star. The pulsation periods of
the high frequencies are in the range of $P_{\rm pul}$ = 0.021$-$0.027 days.
Applying the Table 5 parameters to the relation of $\log Q=\log P_{\rm
pul}+0.5\log g_{1}+0.1M_{\rm bol,1}+\log T_{\rm eff,1}-6.454$ (Petersen &
Jørgensen 1972), we obtained the pulsation constants of $Q$ = 0.014$-$0.018
days. The observed $P_{\rm pul}$ and $Q$ values are typical of $\delta$ Sct
oscillations (Breger 2000; Aerts et al. 2010). Thus, the detached primary
component of X Tri is an oEA pulsator. Our pulsation data corresponding to the
$p$ modes of $\delta$ Sct variables help to probe the stellar envelope of X
Tri A through detailed asteroseismic modeling. However, it is difficult to
probe the deep stellar interiors near the core region and to determine if
there is a phase-dependence in these frequencies, because they did not detect
in $g$-mode pulsations and were observed in only a single TESS band.

Recently, Chen et al. (2022) found that three short-period binaries whose
pulsations were discovered by TESS considerably exceed the upper limit of
$P_{\rm pul}/P_{\rm orb}$ = 0.09$\pm$0.02, meaning that there is no
correlation between them in these stars. X Tri has a relatively shorter
orbital period than other normal oEA binaries (Mkrtichian et al. 2018). Its
pulsation and orbital period ratios are $P_{\rm pul}/P_{\rm orb}<0.008$. The
physical parameters of our target star follow well the empirical relations of
the $\log P_{\rm orb}-\log P_{\rm pul}$, $\log P_{\rm pul}-\log(F/M_{1})$, and
$\log P_{\rm pul}-\log g_{1}$ for the other oEA star with similiar orbital
periods, where $F/M_{1}$ denotes the gravitational pulls onto the pulsating
components exerted by companions (Soydugan et al. 2006; Liakos 2017; Liakos &
Niarchos 2017). Also, our target exhibit almost the same pulsating
characteristics as single $\delta$ Sct pulsators.

The authors would like to thank the BOAO staffs for assistance during our BOES
observations. We also thank the anonymous referee for the careful reading and
valuable comments. This paper includes data collected by the TESS mission and
obtained from MAST. Funding for the TESS mission is provided by the NASA
Explorer Program. IRAF is distributed by the National Optical Astronomy
Observatory, which is operated by the Association of Universities for Research
in Astronomy (AURA) under a cooperative agreement with the National Science
Foundation. This research has made use of the Simbad database maintained at
CDS, Strasbourg, France, and was supported by the KASI grant 2022-1-830-04.
K.H. was supported by the grants 2020R1A4A2002885 and 2022R1I1A1A01053320 from
the National Research Foundation (NRF) of Korea. The research of M.W. was
partially supported by the project Cooperatio - Physics of the Charles
University in Prague.

## References

  * Aerts et al (2010) Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology (Dordrecht: Springer)
  * Bohlin et al (2017) Bohlin, R. C., Mészáros, S., Fleming, S. W., et al. 2017, AJ, 153, 234
  * Bowman et al (2019) Bowman, D. M., Johnston, C., Tkachenko, A., et al. 2019, ApJL, 883, L26
  * Bozkurt et al (1976) Bozkurt, Ş., İbanoǧlu, C., Gülmen, Ö., & Güdür, N. 1976, A&AS, 23, 439
  * Breger (2000) Breger, M. 2000, in ASP Conf. Ser. 210, Delta Scuti and Related Stars, ed. M. Breger and M. Montgomery (San Francisco: ASP), 3
  * Breger et al (1993) Breger, M., Stich, J., Garrido, R., et al. 1993, A&A, 271, 482
  * Chen et al (2022) Chen, X., Ding, X., Cheng, L., et al. 2022, ApJS, 263, 34
  * Duchene & Kraus (2013) Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269
  * Flower (1996) Flower, P. J. 1996, ApJ, 469, 355
  * Frieboes-Conde & Herczeg (1973) Frieboes-Conde, H., & Herczeg, T. 1973, A&AS, 12, 1
  * Gadomski (1932) Gadomski, J. 1932, Warsaw Univ. Reprint No. 15
  * GAIA (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Høg (2000) Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, 27
  * Hong et al (2015) Hong, K., Lee, J. W., Kim, S.-L., et al. 2015, AJ, 150, 131
  * Hong et al (2017) Hong, K., Lee, J. W., Koo, J.-R., et al. 2017, AJ, 153, 247
  * Ibanoǧlu et al (2006) İbanoǧlu, C., Soydugan, F., Soydugan, E., & Dervişoǧlu, A. 2006, MNRAS, 373, 435
  * Ilijic et al (2004) Ilijić, S., Hensberge, H., Pavlovski, K., & Freyhammer, L. M. 2004, in ASP Conf. Ser. 318, Spectroscopically and Spatially Resolving the Components of the Close Binary Stars, ed. R. Hilditch, H. Hensberge, & K. Pavlovski (San Francisco: ASP), 111
  * Kahraman Alicavus et al (2022) Kahraman Aliçavuş, F., Gümüş, D., Kirmizitaş, Ö., et al. 2022, RAA, 22, 085003
  * Kahraman Alicavus et al (2017) Kahraman Aliçavuş, F., Soydugan, E., Smalley, B., & Kubát, J. 2017, MNRAS, 470, 915
  * Kallinger et al (2008) Kallinger, T., Reegen, P., & Weiss, W. W. 2008, A&A, 481, 571
  * Kallrath (2022) Kallrath, J. 2022, Galaxies, 10, 17
  * Kim et al (2022) Kim, H.-Y., Hong, K., Kim, C.-H., et al. 2022, AJ, in press (arXiv:2209.09041)
  * Kim et al (2007) Kim, K.-M., Han, I., Valyavin, G. G., et al. 2007, PASP, 119, 1052
  * Kim et al (2003) Kim, S.-L., Lee, J. W., Kwon, S.-G., Youn, J.-H., Mkrtichian, D. E., & Kim, C. 2003, A&A, 405, 231
  * Kwee & van Woerden (1956) Kwee, K. K., & van Woerden, H. 1956, BAN, 12, 327
  * Lampens (2021) Lampens, P. 2021, Galaxies, 9, 28
  * Lee et al (2021) Lee, J. W., Hong, K., & Kim, H.-Y. 2021, AJ, 161, 129
  * Lee et al (2022a) Lee, J. W., Hong, K., Kim, H.-Y., & Park, J.-H. 2022a, MNRAS, 515, 4702
  * Lee et al (2016) Lee, J. W., Hong, K., Kim, S.-L., & Koo, J.-R. 2016, MNRAS, 460, 4220
  * Lee et al (2018) Lee, J. W., Hong, K., Koo, J.-R., & Park, J.-H. 2018, AJ, 155, 5
  * Lee et al (2020a) Lee, J. W., Hong, K., Koo, J.-R., & Park, J.-H. 2020a, AJ, 159, 24
  * Lee et al (2022b) Lee, J. W., Hong, K., & Park, J.-H. 2022b, MNRAS, 511, 654
  * Lee et al (2014) Lee, J. W., Kim, S.-L., Hong, K., Lee, C.-U., & Koo, J.-R. 2014, AJ, 148, 37
  * Lee et al (2020b) Lee, J. W., Koo, J.-R., Hong, K., & & Park, J.-H. 2020b, AJ, 160, 49
  * Lee et al (2019) Lee, J. W., Kristiansen, M., & Hong, K. 2019, AJ, 157, 223
  * Lenz & Breger (2005) Lenz, P., & Breger, M. 2005, Comm. Asteroseismology, 146, 53
  * Liakos (2017) Liakos, A. 2017, A&A, 607, A85
  * Liakos & Niarchos (2009) Liakos, A., & Niarchos, P. 2009, CoAst, 160, 2
  * Liakos & Niarchos (2015) Liakos, A., & Niarchos, P. 2015, in ASP Conf. Ser. 496, Living Together: Planets, Host Stars and Binaries, ed. S. M. Rucinski, G. Torres, & M. Zejda (San Francisco: ASP), 195
  * Liakos & Niarchos (2017) Liakos, A., & Niarchos, P. 2017, MNRAS, 465, 1181
  * Liakos et al (2010) Liakos, A., Zasche, P., & Niarchos, P. 2010, in ASP Conf. Ser. 435, Binaries - Key to Comprehension of the Universe, ed. A. Prša, & M. Zejda (San Francisco: ASP), 101
  * Mezzetti et al (1980) Mezzetti, M., Cester, B., Giuricin, G., & Mardirossian, F. 1980, A&AS, 39, 265
  * Mkrtichian et al (2002) Mkrtichian, D. E., Kusakin, A. V., Gamarova, A. Y., et al. 2002, in ASP Conf. Ser. 256, Observational Aspects of Pulsating B- and A Stars, ed. C. Sterken, & D. W. Kurtz (San Francisco: ASP), 259
  * Mkrtichian et al (2004) Mkrtichian, D. E., Kusakin, A. V., Rodriguez, E., et al. 2004, A&A, 419, 1015
  * Mkrtichian et al (2018) Mkrtichian, D. E., Lehmann, H., Rodríguez, E., et al. 2018, MNRAS, 475, 4745
  * Murphyz (2018) Murphy, S. J. 2018, preprint (arXiv:1811.12659)
  * Neujmin (1922) Neujmin, G. 1922, AN, 217, 192
  * Pecaut & Mamajek (2013) Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9
  * Petersen & Jorgensen (1972) Petersen, J. O., & Jørgensen, H. E. 1972, A&A, 17, 367
  * Qian (2002) Qian, S. 2002, PASP, 114, 650
  * Rafert (1982) Rafert, J. B. 1982, PASP, 94, 485
  * Ricker et al (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, JATIS, 1, 014003
  * Rovithis-Livaniou et al (2000) Rovithis-Livaniou, H., Kranidiotis, A. N., Rovithis, P., & Athanassiades, G. 2000, A&A, 354, 904
  * Shi et al (2022) Shi, X.-D., Qian, S.-B., & Li, L.-J. 2022, ApJS, 259, 50
  * Soydugan et al (2006) Soydugan, E., İbanoǧlu, C., Soydugan, F., Akan, M. C., & Demircan, O. 2006, MNRAS, 366, 1298
  * Struve (1946) Struve, O. 1946, ApJ, 104, 253
  * Torres (2010) Torres, G. 2010, AJ, 140, 1158
  * Tout & Hall (1991) Tout, C. A., & Hall, D. S. 1991, MNRAS, 253, 9
  * Turner et al (2014) Turner, G., Katichuck, R., & Holaday, J. 2014, JAVSO, 42, 134
  * Van Hamme (1993) Van Hamme, W. 1993, AJ, 106, 209
  * Van Hamme & Wilson (1990) Van Hamme, W., & Wilson, R. E. 1990, AJ, 100, 198
  * Van Hamme & Wilson (2007) Van Hamme, W., & Wilson, R. E. 2007, ApJ, 661, 1129
  * Walker (1921) Walker, A. D. 1921, HarCi, 225, 1
  * Wilson & Devinney (1971) Wilson, R. E., & Devinney, E. J. 1971, ApJ, 166, 605
  * Wood (1950) Wood, F. B. 1950, ApJ, 112, 196
  * Wood & Forbes (1963) Wood, B. D., & Forbes, J. E. 1963, AJ, 68, 257
  * Zhang et al (2013) Zhang, X. B., Luo, C. Q., & Fu, J. N. 2013, ApJ, 777, 77
  * Zhou (2010) Zhou, A. Y. 2010, arXiv:1002.2729v6

Figure 1: Light curve of X Tri with the fitted model. In the top panel, the
gray circles are individual measurements from TESS, and the red solid curve is
computed with a cool spot on the secondary star. The light curve residuals
corresponding to the unspotted and cool-spot models are plotted in the second
and third panels, respectively. The normal points phase-binning at intervals
of 0.001 are plotted as the cyan dots in the top panel and displaced
vertically for clarity. The residuals from the mean curve are shown in the
bottom panel.  Figure 2: Historical $O-C$ diagram of X Tri since discovery,
covering over 100 years of observations. The individual primary minima are
denoted by blue circles and the secondary by yellow triangles. Large symbols
correspond to more precise photoelectric or CCD measurements. Our new TESS
minima are plotted as a bulk of points near the epoch 11800. In the upper
panel constructed with the linear terms of equation (2), the quadratic fit is
given by the red parabolic curve. The timing residuals after subtracting the
parabolic trend are plotted in the lower panel. The quasi-sinusoidal
variations of $O-C$ values up to 0.01 days are clearly visible.  Figure 3: The
upper panel displays the trailed spectra of X Tri in Ca I doublet 6103/6122 Å.
In the lower panel, the gray circles and blue and red lines represent the
observed spectrum at an orbital phase of 0.73 and the synthetic spectra of the
primary (X Tri A; $T_{\rm eff,1}$ = 7900 K, $\log$ $g_{1}$ = 4.3,
$v_{1}$$\sin$$i$ = 84 km s-1) and secondary (X Tri B; $T_{\rm eff,2}$ = 4990
K, $\log$ $g_{2}$ = 3.9, $v_{2}$$\sin$$i$ = 103 km s-1) components,
respectively.  Figure 4: RV curves of X Tri with fitted models. The blue and
green symbols are the primary and secondary measurements, respectively. The
red solid curves represent the results from a consistent light and RV curve
analysis with the W-D code. The gray line in the upper panel represents the
system velocity of $-$9.21 km s-1. The lower panel shows the RV residuals
between observations and models.  Figure 5: Seven spectral regions of the
primary star. The gray circles represent the disentangling spectrum obtained
by the FDBinary code. The red lines denote the synthetic spectra of our best-
fit parameters $T_{\rm eff,1}$ = 7900 K and $v_{1}$$\sin$$i$ = 84 km s-1.
Figure 6: Variations in light levels for X Tri at four characteristic phases:
Max I, Min I, Max II and Min II. See text in detail.  Figure 7: Light curve
residuals from 1000 normal points distributed in BJD. The vertical dotted
lines indicate the primary minimum epochs. The lower panel presents a short
section of the residuals marked using the inset box in the upper panel. The
synthetic curve is computed from the 16-frequency fit to the out-of-eclipse
part of the residuals.  Figure 8: Amplitude spectra before (top panel) and
after pre-whitening of the first 4 frequencies (middle panel) and all 16
frequencies (bottom panel) from the PERIOD04 program. The red line in the
bottom panel corresponds to four times the noise spectrum.

Table 1: TESS Eclipse Timings for X Tri  BJD | Error | Min | BJD | Error | Min
---|---|---|---|---|---
2,458,765.14863 |  $\pm$0.00001 | I | 2,458,779.23707 |  $\pm$0.00010 | II
2,458,765.63585 |  $\pm$0.00011 | II | 2,458,779.72151 |  $\pm$0.00001 | I
2,458,766.12013 |  $\pm$0.00001 | I | 2,458,780.20852 |  $\pm$0.00010 | II
2,458,766.60747 |  $\pm$0.00013 | II | 2,458,780.69301 |  $\pm$0.00001 | I
2,458,767.09168 |  $\pm$0.00001 | I | 2,458,781.18010 |  $\pm$0.00009 | II
2,458,767.57895 |  $\pm$0.00014 | II | 2,458,781.66456 |  $\pm$0.00001 | I
2,458,768.06322 |  $\pm$0.00001 | I | 2,458,782.15148 |  $\pm$0.00010 | II
2,458,768.55024 |  $\pm$0.00009 | II | 2,458,782.63613 |  $\pm$0.00001 | I
2,458,769.03474 |  $\pm$0.00001 | I | 2,458,783.12299 |  $\pm$0.00009 | II
2,458,769.52184 |  $\pm$0.00012 | II | 2,458,783.60761 |  $\pm$0.00001 | I
2,458,770.00627 |  $\pm$0.00001 | I | 2,458,784.09455 |  $\pm$0.00010 | II
2,458,770.49336 |  $\pm$0.00011 | II | 2,458,784.57913 |  $\pm$0.00001 | I
2,458,770.97779 |  $\pm$0.00001 | I | 2,458,785.06601 |  $\pm$0.00010 | II
2,458,771.46491 |  $\pm$0.00011 | II | 2,458,785.55067 |  $\pm$0.00001 | I
2,458,771.94933 |  $\pm$0.00001 | I | 2,458,786.03754 |  $\pm$0.00009 | II
2,458,772.43605 |  $\pm$0.00005 | II | 2,458,786.52221 |  $\pm$0.00001 | I
2,458,777.77844 |  $\pm$0.00001 | I | 2,458,787.00902 |  $\pm$0.00008 | II
2,458,778.26567 |  $\pm$0.00013 | II | 2,458,787.49374 |  $\pm$0.00001 | I
2,458,778.74997 |  $\pm$0.00001 | I | 2,458,787.98039 |  $\pm$0.00007 | II
Table 2: Radial velocities of X Tri  |  | Fe I $\lambda$4957 |  | Ca I $\lambda$6103
---|---|---|---|---
BJD | UT date | $V_{1}$ | $\sigma_{1}$ | $V_{2}$ | $\sigma_{2}$ |  | $V_{1}$ | $\sigma_{1}$ | $V_{2}$ | $\sigma_{2}$
(2,457,000+) |  | (km s-1) | (km s-1) | (km s-1) | (km s-1) |  | (km s-1) | (km s-1) | (km s-1) | (km s-1)
1,820.0618 | 2019-12-02 | $0.1$ | 6.1 |  |  |  | $4.8$ | 2.2 |  |
1,820.0816 | 2019-12-02 | $12.1$ | 4.4 |  |  |  | $19.9$ | 6.2 |  |
1,820.1013 | 2019-12-02 | $27.9$ | 4.9 |  |  |  | $31.6$ | 6.3 |  |
1,820.1210 | 2019-12-02 | $41.4$ | 5.9 |  |  |  | $41.8$ | 6.7 |  |
1,820.1408 | 2019-12-02 | $54.4$ | 7.6 |  |  |  | $53.6$ | 5.0 | $-136.8$ | 8.5
1,820.1605 | 2019-12-02 | $64.3$ | 5.6 | $-164.5$ | 10.0 |  | $64.3$ | 7.7 | $-145.2$ | 10.7
1,820.1803 | 2019-12-02 | $71.3$ | 4.5 | $-176.0$ | 8.2 |  | $73.8$ | 4.1 | $-171.9$ | 12.4
2,145.1004 | 2020-10-22 | $-72.8$ | 8.8 |  |  |  | $-70.2$ | 6.7 |  |
2,145.1215 | 2020-10-22 | $-73.4$ | 7.4 |  |  |  | $-88.0$ | 11.9 |  |
2,145.2379 | 2020-10-22 | $-115.3$ | 6.4 | $192.6$ | 9.0 |  | $-116.4$ | 7.0 | $195.5$ | 6.7
2,145.2590 | 2020-10-22 | $-120.6$ | 3.1 | $195.6$ | 7.0 |  | $-122.8$ | 7.3 | $189.7$ | 9.0
2,145.2801 | 2020-10-22 | $-116.9$ | 4.5 | $193.7$ | 4.5 |  | $-111.8$ | 5.3 | $193.5$ | 5.5
2,145.3013 | 2020-10-22 | $-114.1$ | 5.4 | $191.5$ | 5.5 |  | $-122.9$ | 6.8 | $193.5$ | 9.1
2,145.3224 | 2020-10-22 | $-110.5$ | 6.8 | $180.4$ | 8.2 |  | $-106.0$ | 4.3 | $181.2$ | 5.3
2,146.1355 | 2020-10-23 | $-98.6$ | 3.1 | $166.4$ | 6.1 |  | $-100.2$ | 9.7 | $143.9$ | 13.5
2,146.1567 | 2020-10-23 | $-107.8$ | 6.4 | $173.3$ | 6.8 |  | $-107.2$ | 5.4 | $163.4$ | 10.2
2,146.1778 | 2020-10-23 | $-112.8$ | 3.9 | $183.0$ | 4.1 |  | $-111.3$ | 5.2 | $177.1$ | 10.3
2,146.1989 | 2020-10-23 | $-117.5$ | 5.4 | $195.5$ | 9.6 |  | $-117.5$ | 4.5 | $187.6$ | 8.5
2,146.2201 | 2020-10-23 | $-120.6$ | 5.7 | $191.1$ | 8.4 |  | $-121.6$ | 6.3 | $186.7$ | 9.6
2,146.2412 | 2020-10-23 | $-117.7$ | 5.0 | $197.1$ | 5.5 |  | $-117.1$ | 4.2 | $190.5$ | 7.2
2,146.2623 | 2020-10-23 | $-116.7$ | 2.8 | $195.3$ | 4.5 |  | $-118.9$ | 4.0 | $184.1$ | 10.0
2,146.2835 | 2020-10-23 | $-112.2$ | 5.0 | $191.8$ | 5.1 |  | $-115.9$ | 5.2 | $183.0$ | 7.6
2,146.3046 | 2020-10-23 | $-105.5$ | 4.0 | $180.0$ | 5.9 |  | $-102.3$ | 9.7 | $169.4$ | 17.6
2,146.3257 | 2020-10-23 | $-101.3$ | 2.6 |  |  |  | $-101.0$ | 5.6 | $168.9$ | 7.8
2,167.0020 | 2020-11-13 | $67.7$ | 3.4 |  |  |  | $74.3$ | 5.6 | $-147.0$ | 16.9
2,167.0230 | 2020-11-13 | $75.3$ | 3.0 | $-167.8$ | 9.8 |  | $77.9$ | 5.2 | $-175.2$ | 17.1
2,167.0441 | 2020-11-13 | $84.1$ | 5.2 | $-179.7$ | 5.0 |  | $85.4$ | 3.7 | $-191.5$ | 9.8
2,167.0654 | 2020-11-13 | $87.6$ | 5.4 | $-204.1$ | 6.6 |  | $91.1$ | 6.0 | $-190.8$ | 9.1
2,167.0868 | 2020-11-13 | $94.4$ | 5.5 | $-197.5$ | 9.0 |  | $92.0$ | 4.9 | $-201.7$ | 12.3
2,167.1078 | 2020-11-13 | $96.7$ | 4.5 | $-207.7$ | 7.3 |  | $96.2$ | 6.7 | $-200.8$ | 8.7
2,169.0452 | 2020-11-15 | $94.7$ | 4.2 | $-205.1$ | 8.0 |  | $90.4$ | 6.0 | $-207.0$ | 22.0
2,169.0664 | 2020-11-15 | $94.2$ | 5.0 | $-212.0$ | 7.8 |  | $98.6$ | 6.5 | $-205.4$ | 9.6
2,169.0883 | 2020-11-15 | $93.3$ | 5.7 | $-204.4$ | 6.7 |  | $98.0$ | 5.7 | $-209.4$ | 8.0
2,169.1095 | 2020-11-15 | $91.3$ | 3.0 | $-200.2$ | 8.2 |  | $92.8$ | 7.1 | $-206.1$ | 8.2
2,169.1308 | 2020-11-15 | $85.0$ | 4.4 | $-190.9$ | 7.7 |  | $87.2$ | 5.3 | $-193.4$ | 10.1
Table 3: Light levels of X Tri at four different phases  Mean Time | Min I | Max I | Min II | Max II
---|---|---|---|---
(BJD) | (mag) | (mag) | (mag) | (mag)
2,458,765.17500 | 10.0972$\pm$0.0252 | 8.5652$\pm$0.0013 | 8.8562$\pm$0.0045 | 8.5708$\pm$0.0014
2,458,766.14726 | 10.0985$\pm$0.0275 | 8.5660$\pm$0.0014 | 8.8564$\pm$0.0040 | 8.5727$\pm$0.0020
2,458,767.11951 | 10.0986$\pm$0.0258 | 8.5657$\pm$0.0019 | 8.8568$\pm$0.0046 | 8.5713$\pm$0.0014
2,458,768.09176 | 10.0974$\pm$0.0263 | 8.5654$\pm$0.0024 | 8.8567$\pm$0.0039 | 8.5713$\pm$0.0013
2,458,769.06401 | 10.0986$\pm$0.0254 | 8.5650$\pm$0.0021 | 8.8554$\pm$0.0033 | 8.5725$\pm$0.0017
2,458,770.03626 | 10.1001$\pm$0.0268 | 8.5666$\pm$0.0011 | 8.8562$\pm$0.0036 | 8.5727$\pm$0.0011
2,458,771.00851 | 10.0983$\pm$0.0252 | 8.5656$\pm$0.0024 | 8.8557$\pm$0.0040 | 8.5720$\pm$0.0008
2,458,771.98006 | 10.0955$\pm$0.0265 | 8.5657$\pm$0.0019 | 8.8565$\pm$0.0034 | 8.5722$\pm$0.0008
2,458,778.22116 | 10.0974$\pm$0.0255 | 8.5660$\pm$0.0013 | 8.8562$\pm$0.0044 | 8.5719$\pm$0.0008
2,458,779.19340 | 10.0964$\pm$0.0254 | 8.5658$\pm$0.0006 | 8.8569$\pm$0.0042 | 8.5715$\pm$0.0016
2,458,780.16563 | 10.0962$\pm$0.0256 | 8.5660$\pm$0.0013 | 8.8564$\pm$0.0037 | 8.5705$\pm$0.0013
2,458,781.13787 | 10.0957$\pm$0.0262 | 8.5654$\pm$0.0010 | 8.8555$\pm$0.0026 | 8.5716$\pm$0.0017
2,458,782.11010 | 10.0962$\pm$0.0248 | 8.5654$\pm$0.0012 | 8.8561$\pm$0.0037 | 8.5718$\pm$0.0011
2,458,783.08232 | 10.0916$\pm$0.0260 | 8.5653$\pm$0.0012 | 8.8578$\pm$0.0033 | 8.5712$\pm$0.0010
2,458,784.05455 | 10.0931$\pm$0.0254 | 8.5648$\pm$0.0019 | 8.8565$\pm$0.0035 | 8.5717$\pm$0.0005
2,458,785.02678 | 10.0931$\pm$0.0257 | 8.5647$\pm$0.0010 | 8.8553$\pm$0.0025 | 8.5709$\pm$0.0008
2,458,785.99900 | 10.0926$\pm$0.0249 | 8.5651$\pm$0.0013 | 8.8567$\pm$0.0043 | 8.5721$\pm$0.0012
2,458,786.97122 | 10.0920$\pm$0.0258 | 8.5648$\pm$0.0007 | 8.8565$\pm$0.0042 | 8.5719$\pm$0.0009
2,458,787.90941 | 10.0879$\pm$0.0254 | 8.5660$\pm$0.0015 | 8.8564$\pm$0.0043 | 8.5709$\pm$0.0011
Average | 10.0956$\pm$0.0031 | 8.5655$\pm$0.0005 | 8.8563$\pm$0.0006 | 8.5717$\pm$0.0006
Table 4: Light and RV Parameters of X Tri

Parameter | Primary | Secondary
---|---|---
$T_{0}$ (BJD) | 2,458,765.148597$\pm$0.000028
$P_{\rm orb}$ (day) | 0.9715319$\pm$0.0000018
$a$ ($R_{\odot}$) | 6.105$\pm$0.021
$\gamma$ (km s-1) |  $-$9.21$\pm$0.29
$K_{1}$ (km s-1) | 108.15$\pm$0.50
$K_{2}$ (km s-1) | 209.83$\pm$0.98
$q$ | 0.5154$\pm$0.0035
$i$ (deg) | 88.90$\pm$0.14
$T_{\rm eff}$ (K) | 7900$\pm$110 | 4990$\pm$50
$\Omega$ | 4.220$\pm$0.014 | 2.905
$X_{\rm bol}$, $Y_{\rm bol}$ | 0.673, 0.203 | 0.642, 0.165
$x_{T_{\rm P}}$, $y_{T_{\rm P}}$ | 0.574$\pm$0.031, 0.214 | 0.695$\pm$0.023, 0.188
$l$/($l_{1}$+$l_{2}$+$l_{3}$) | 0.7210$\pm$0.0031 | 0.2416
$l_{3}$a | 0.0374$\pm$0.0024
$r$ (pole) | 0.2687$\pm$0.0013 | 0.3022$\pm$0.0011
$r$ (point) | 0.2786$\pm$0.0015 | 0.4323$\pm$0.0013
$r$ (side) | 0.2728$\pm$0.0013 | 0.3154$\pm$0.0012
$r$ (back) | 0.2766$\pm$0.0014 | 0.3479$\pm$0.0012
$r$ (volume)b | 0.2728$\pm$0.0014 | 0.3232$\pm$0.0012
Spot Parameters:
Colatitude (deg) |  | 59.9$\pm$2.4
Longitude (deg) |  | 335.9$\pm$1.9
Radius (deg) |  | 16.8$\pm$2.8
$T$spot/$T$local |  | 0.885$\pm$0.039

aafootnotetext: Value at 0.25 orbital phase.

bbfootnotetext: Mean volume radius.

Table 5: Absolute Parameters of X Tri  Parameter | Primary | Secondary
---|---|---
$M$ ($M_{\odot}$) | 2.137$\pm$0.018 | 1.101$\pm$0.010
$R$ ($R_{\odot}$) | 1.664$\pm$0.010 | 1.972$\pm$0.010
$\log$ $g$ (cgs) | 4.325$\pm$0.006 | 3.890$\pm$0.006
$\rho$ ($\rho_{\odot}$) | 0.464$\pm$0.010 | 0.1439$\pm$0.002
$v_{\rm sync}$ (km s-1) | 86.71$\pm$0.54 | 102.72$\pm$0.52
$v$$\sin$$i$ (km s-1) | 84$\pm$6 |
$T_{\rm eff}$ (K) | 7900$\pm$110 | 4990$\pm$50
$L$ ($L_{\odot}$) | 9.67$\pm$0.55 | 2.16$\pm$0.09
$M_{\rm bol}$ (mag) | 2.267$\pm$0.062 | 3.894$\pm$0.045
BC (mag) | 0.022$\pm$0.001 |  $-$0.292$\pm$0.017
$M_{\rm V}$ (mag) | 2.245$\pm$0.062 | 4.186$\pm$0.048
Distance (pc) | 215$\pm$6
Table 6: Results of the multiple frequency analysis for X Tria

| Frequency | Amplitude | Phase | SNRb | Remark
---|---|---|---|---|---
| (day-1) | (mmag) | (rad) |  |
$f_{1}$ | 40.6840$\pm$0.0016 | 0.322$\pm$0.062 | 1.55$\pm$0.57 | 8.86 |
$f_{2}$ | 45.2287$\pm$0.0010 | 0.486$\pm$0.057 | 0.14$\pm$0.35 | 14.49 |
$f_{3}$ | 36.1139$\pm$0.0014 | 0.320$\pm$0.054 | 1.15$\pm$0.49 | 10.16 | $2f_{1}-f_{2}$
$f_{4}$ | 43.5160$\pm$0.0014 | 0.352$\pm$0.059 | 0.06$\pm$0.49 | 10.19 |
$f_{5}$ | 0.6402$\pm$0.0032 | 0.253$\pm$0.095 | 5.56$\pm$1.10 | 4.57 |
$f_{6}$ | 37.3009$\pm$0.0017 | 0.269$\pm$0.055 | 2.28$\pm$0.60 | 8.39 |
$f_{7}$ | 43.8975$\pm$0.0019 | 0.266$\pm$0.059 | 2.10$\pm$0.64 | 7.78 |
$f_{8}$ | 45.2902$\pm$0.0015 | 0.330$\pm$0.058 | 4.25$\pm$0.51 | 9.82 | $2f_{1}-f_{3}$
$f_{9}$ | 38.6236$\pm$0.0015 | 0.321$\pm$0.056 | 2.74$\pm$0.51 | 9.76 |
$f_{10}$ | 51.8274$\pm$0.0017 | 0.212$\pm$0.042 | 1.13$\pm$0.58 | 8.60 | $2f_{2}-f_{9}$
$f_{11}$ | 42.2950$\pm$0.0025 | 0.209$\pm$0.063 | 1.11$\pm$0.88 | 5.72 |
$f_{12}$ | 0.9539$\pm$0.0041 | 0.193$\pm$0.094 | 4.57$\pm$1.43 | 3.50 |
$f_{13}$ | 47.2891$\pm$0.0015 | 0.298$\pm$0.054 | 5.47$\pm$0.53 | 9.43 |
$f_{14}$ | 2.0964$\pm$0.0036 | 0.207$\pm$0.088 | 1.62$\pm$1.25 | 4.01 | $2f_{\rm orb}$
$f_{15}$ | 0.5066$\pm$0.0040 | 0.198$\pm$0.094 | 6.06$\pm$1.39 | 3.61 |
$f_{16}$ | 39.7789$\pm$0.0034 | 0.147$\pm$0.059 | 1.17$\pm$1.18 | 4.24 | $2f_{4}-f_{13}$

aafootnotetext: Frequencies are listed in order of detection.

bbfootnotetext: Calculated in a range of 5 day-1 around each frequency.