A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography <s> Optimization-based image reconstruction <s> List of Algorithms Preface 1. Introduction. Brief Overview of the State of the Art Notation Review of Relevant Linear Algebra Part I. Krylov Subspace Approximations. 2. Some Iteration Methods. Simple Iteration Orthomin(1) and Steepest Descent Orthomin(2) and CG Orthodir, MINRES, and GMRES Derivation of MINRES and CG from the Lanczos Algorithm 3. Error Bounds for CG, MINRES, and GMRES. Hermitian Problems-CG and MINRES Non-Hermitian Problems-GMRES 4. Effects of Finite Precision Arithmetic. Some Numerical Examples The Lanczos Algorithm A Hypothetical MINRES/CG Implementation A Matrix Completion Problem Orthogonal Polynomials 5. BiCG and Related Methods. The Two-Sided Lanczos Algorithm The Biconjugate Gradient Algorithm The Quasi-Minimal Residual Algorithm Relation Between BiCG and QMR The Conjugate Gradient Squared Algorithm The BiCGSTAB Algorithm Which Method Should I Use? 6. Is There A Short Recurrence for a Near-Optimal Approximation? The Faber and Manteuffel Result Implications 7. Miscellaneous Issues. Symmetrizing the Problem Error Estimation and Stopping Criteria Attainable Accuracy Multiple Right-Hand Sides and Block Methods Computer Implementation Part II. Preconditioners. 8. Overview and Preconditioned Algorithms. 9. Two Example Problems. The Diffusion Equation The Transport Equation 10. Comparison of Preconditioners. Jacobi, Gauss--Seidel, SOR The Perron--Frobenius Theorem Comparison of Regular Splittings Regular Splittings Used with the CG Algorithm Optimal Diagonal and Block Diagonal Preconditioners 11. Incomplete Decompositions. Incomplete Cholesky Decomposition Modified Incomplete Cholesky Decomposition 12. Multigrid and Domain Decomposition Methods. Multigrid Methods Basic Ideas of Domain Decomposition Methods. <s> BIB001 </s> A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography <s> Optimization-based image reconstruction <s> PET and SPECT are two of todays most important medical-imaging methods, providing images that reveal subtle information about physiological processes in humans and animals. "Emission Tomography: The Fundamentals of PET and SPECT" explains the physics and engineering principles of these important functional-imaging methods. The technology of emission tomography is covered in detail, including historical origins, scientific and mathematical foundations, imaging systems and their components, image reconstruction and analysis, simulation techniques, and clinical and laboratory applications. The book describes the state of the art of emission tomography, including all facets of conventional SPECT and PET, as well as contemporary topics, such as iterative image reconstruction, small-animal imaging, and PET/CT systems. This book is intended as a textbook and reference resource for graduate students, researchers, medical physicists, biomedical engineers, and professional engineers and physicists in the medical-imaging industry. Thorough tutorials of fundamental and advanced topics are presented by dozens of the leading researchers in PET and SPECT. SPECT has long been a mainstay of clinical imaging, and PET is now one of the worlds fastest growing medical imaging techniques, owing to its dramatic contributions to cancer imaging and other applications. "Emission Tomography: The Fundamentals of PET and SPECT" is an essential resource for understanding the technology of SPECT and PET, the most widely used forms of molecular imaging. The features are: contains thorough tutorial treatments, coupled with coverage of advanced topics; three of the four holders of the prestigious Institute of Electrical and Electronics Engineers Medical Imaging Scientist Award are chapter contributors; and includes color artwork. <s> BIB002 </s> A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography <s> Optimization-based image reconstruction <s> In this paper we analyze several new methods for solving optimization problems with the objective function formed as a sum of two convex terms: one is smooth and given by a black-box oracle, and another is general but simple and its structure is known. Despite to the bad properties of the sum, such problems, both in convex and nonconvex cases, can be solved with efficiency typical for the good part of the objective. For convex problems of the above structure, we consider primal and dual variants of the gradient method (converge as O (1/k)), and an accelerated multistep version with convergence rate O (1/k2), where k isthe iteration counter. For all methods, we suggest some efficient "line search" procedures and show that the additional computational work necessary for estimating the unknown problem class parameters can only multiply the complexity of each iteration by a small constant factor. We present also the results of preliminary computational experiments, which confirm the superiority of the accelerated scheme. <s> BIB003 </s> A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography <s> Optimization-based image reconstruction <s> This book presents the state of the art in sparse and multiscale image and signal processing, covering linear multiscale transforms, such as wavelet, ridgelet, or curvelet transforms, and non-linear multiscale transforms based on the median and mathematical morphology operators. Recent concepts of sparsity and morphological diversity are described and exploited for various problems such as denoising, inverse problem regularization, sparse signal decomposition, blind source separation, and compressed sensing. This book weds theory and practice in examining applications in areas such as astronomy, biology, physics, digital media, and forensics. A final chapter explores a paradigm shift in signal processing, showing that previous limits to information sampling and extraction can be overcome in very significant ways. Matlab and IDL code accompany these methods and applications to reproduce the experiments and illustrate the reasoning and methodology of the research available for download at the associated Web site. <s> BIB004 </s> A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography <s> Optimization-based image reconstruction <s> "Fixed-Point Algorithms for Inverse Problems in Science and Engineering" presents some ofthe most recent work from top-notch researchers studying projection and other first-order fixed-point algorithms in several areas of mathematics and the applied sciences. The material presented provides a survey of the state-of-the-art theory and practice in fixed-point algorithms, identifying emerging problems driven by applications, and discussing new approaches for solving these problems. This book incorporates diverse perspectives from broad-ranging areas of research including, variational analysis, numerical linear algebra, biotechnology, materials science, computational solid-state physics, and chemistry. Topics presented include: Theory of Fixed-point algorithms: convex analysis, convex optimization, subdifferential calculus, nonsmooth analysis, proximal point methods, projection methods, resolvent and related fixed-point theoretic methods, and monotone operator theory. Numerical analysis of fixed-point algorithms: choice of step lengths, of weights, of blocks for block-iterative and parallel methods, and of relaxation parameters; regularization of ill-posed problems; numerical comparison of various methods. Areas of Applications: engineering (image and signal reconstruction and decompression problems), computer tomography and radiation treatment planning (convex feasibility problems), astronomy (adaptive optics), crystallography (molecular structure reconstruction), computational chemistry (molecular structure simulation) and other areas. Because of the variety of applications presented, this book can easily serve as a basis for new and innovated research and collaboration. <s> BIB005
A general approach to PACT image reconstruction is to formulate the sought-after estimate of p 0 (r) as the solution of an optimization problem. In fact, most modern image reconstruction methods for computed imaging modalities including X-ray computed tomography and magnetic resonance imaging are formulated in this way (Fessler, 1994; BIB002 Pan et al., 2009 ). There are potential practical and conceptual advantages to optimization-based image reconstruction over analytic methods. For example, because they are based on D-D forward models, optimizationbased image reconstruction methods can comprehensively compensate for the imaging physics and other physical factors such as responses of the measurement system that are not easily incorporated into an analytic method. Moreover, such methods provide a general framework for incorporating regularization, which can mitigate the effects of measurement noise and data incompleteness. Because optimization-based methods are often implemented by use of iterative algorithms, they are generally more computationally demanding than analytic methods; however the use of modern parallel computing technologies can render iterative three dimensional model-based reconstruction scheme for arbitrary photoacoustic acquisition geometries feasible for many PACT applications. Consider a D-D imaging model u = Hθ as described in Section 3. An optimizationbased image reconstruction method seeks to determine an estimate of θ by solving an optimization program that can be generally specified aŝ Here, F (·) is the objective function to be minimized that depends on the D-D forward operator H, the measured data u, and the vector θ that specifies the estimate of p 0 (r). The functions {f i (θ)} Nc i=0 and the closed set {C i } Nc i=0 describe the set of N c constraints that the solution must satisfy. It is important to note that numerous choices must be made in the design of the reconstruction method. These include the specification of H, F (·), and f (θ), as well as the optimization algorithm that is employed to solve Eqn. (47). It is the joint specification of these quantities that defines the sought after solution and the numerical properties of the reconstruction method (Zhang et al., 2009  Since the initial pressure distribution is non-negative, the above problem can be constrained so the solution satisfies θ ≥ 0. The objective function in the PLS estimator is expressed as two terms. The quantity W is referred to as the data fidelity term that corresponds to a weighted least squares functional. The matrix W ∈ R QK×QK that defines the weighted l 2 norm is symmetric and positive definite, whereas u ∈ R QK×1 corresponds to the measurement vector. This data fidelity functional is convex and differentiable with respect to θ. For the case of Gaussian measurement noise, this data fidelity term can be interpreted as a negative log-likelihood function; however, its use is not restricted to that case. The quantity R(θ) is a penalty term that can be designed to regularize the inverse problem, and γ ∈ R is a regularization parameter that controls the amount of regularization. A classic form of regularization is Tikhonov regularization (Golub et al., 1999) , where the regularization term is specified as where P ∈ R N ×N is a symmetric positive definite matrix. As the Tikhonov regularization term is differentiable with respect to θ, a variety of gradient-based methods can be employed to solve Eqn. (48). Gradient-based methods that can be employed to solve Eqn. (48) can be broadly grouped into two classes: first order methods and second order methods. As part of first order methods, the Taylor approximation used to compute the descent direction involves a linear or first order approximation of the cost function. Hence, information about the gradient of the function is employed to compute the descent direction. Some of the most commonly employed first order methods include the steepest descent method, and the conjugate gradient method. There also exists a class of first order methods that employ information about the past gradients/momentum to speed up the convergence rate of first order algorithms. Such methods are referred as Nesterov methods, whereby linear combinations of present and past gradients are used to compute the descent direction (Nesterov, 1983; Nesterov, 1988) . In addition, a family of methods called Krylov-based methods, the subset of which is are the conjugate gradient methods, are also employed to solve smooth, convex programs defined in Eqn. (48) (Paige and Saunders, 1982; Gutknecht, 2007; BIB001 . The design and study of variants of the Nesterov method and the Krylov-based methods to solve smooth, convex optimization programs is an active research area (Ghai et al., 2019;  Dax, 2019) . The second class of methods regularly employed to optimization programs are the second order methods. The Taylor approximation used to compute the descent direction for second order methods is a quadratic or second order approximation of the cost function. Hence, information about the inverse Hessian of the cost function is employed to compute the descent direction. Methods that explicitly compute the inverse Hessian are referred to as Newton methods. Although the newton methods have favorable converge properties, the computational burden associated with computing a Hessian and inverting it is prohibitively large. Hence, a variety of methods that approximate the inverse Hessian from first order gradient information are commonly employed. Such methods are referred to as Quasi-Newton methods. A variety of Quasi-Newton methods have been proposed and comprehensively studied (Wright and Nocedal, 1999;  Nocedal, 1980; Conn et al., 1991; Khalfan et al., 1993; Dai, 2002) . The set of methods discusssed above can handle cost functions that have smooth regularization and data fidelity terms that are differentiable. However, certain modern approaches to regularization employ non-smooth choices for R(θ) based on the l 1 -norm: where Φ ∈ R N ×N is a sparsifying transformation that is chosen such that θ is sparse in its range. Popular choices for Φ include the discrete wavelet transform or finite difference operators. Total variation (TV) regularization (Sidky and Pan, 2008) can be achieved as a special case when the latter are employed. Such choices are based on the observation that many objects possess a sparse representation in the wavelet or gradient domain BIB004 Chartrand, 2009) . Although Eqn. (50) is convex, it is not differentiable with respect to θ. However, a variety of non-smooth optimization methods, such as proximal-gradient methods, can be employed to solve Eqn. (48) in this case. Proximal methods are a type of forward-backward splitting approach, which alternates between computing a gradient step on the data fidelity term and a solution to a proximal problem involving the regularization (Parikh et al., 2014; BIB005 Beck and Teboulle, 2009b) . Moreover, the proximal gradient methods can also be accelerated by utilizing momentum information from past gradients or by deploying second order information to solve the proximal problem  Beck and Teboulle, 2009b; Beck and Teboulle, 2009a; Pan et al., 2013; BIB003  Stella et al., 2017) . Examples of PACT images of a live mouse reconstructed by use of two different reconstruction methods are displayed in Fig. 3 . A description of the imaging system and experimental details have been described in previous works Brecht et al., 2009 ). The small animal imager consisted of an arc-shaped probe of 64 transducers along the vertical axis. Thus, a tomographic view consisted of data recorded at 64 transducers along the arc. The complete tomographic dataset was acquired by rotating the object 360 • around the vertical axis. A 180-view dataset was used to reconstruct 3D images of the mouse trunk as shown in Fig. 3 . The image shown in Fig. 3(a) was reconstructed by use of a FBP algorithm, while the image in Fig. 3(b) was reconstructed by use of the PLS estimator in Eqn. (48) that employed a TV penalty. The image corresponding to the PLS-TV estimate possesses a higher contrast than the images reconstructed by use of the FBP algorithm and reveal a much sharper body vascular tree. These observations are consistent with the fact that optimization-based image reconstruction methods can produce images that possess different physical and statistical characteristics than the ones produced by analytical method. An important consideration in the formulation of an optimization-based image reconstruction method is the choice of the penalty function. To demonstrate how different choices can influence image quality, Fig. 4 displays images of a simple experimental phantom reconstructed by use of Eqn. (48) with two different choices for the penalty function. Figures 4 (a) and (b) correspond to use of quadratic smoothness Laplacian regularization and TV regularization, respectively. The quadratic Laplacian penalty corresponds to the l 2 -norm of a discrete Laplacian operator acting on the image estimate. In this example, use of the TV penalty resulted in images with reduced artifact levels and enhanced spatial resolution as compared to use of the quadratic Laplacian penalty. All the images shown in Figures 3 and 4 were reconstructed by use of an iterative image reconstruction model that assumed homogeneous acoustic media and thus could not account for the acoustic heterogeneity of the medium. However, knowledge about the spatial distribution of acoustic heterogeneity of the medium can be incorporated into an iterative reconstruction algorithm based on the imaging model described in Eqn. (28). In vivo experimental studies have been conducted to validate the aforementioned iterative algorithm whereby the spatial distribution of the acoustic properties of the medium are (48) with the penalty function specified as the (a) l 2 -norm of a discrete Laplacian and (b) TV semi-norm. Both the images are of the same size and are normalized to the same grayscale window.The images were reproduced from the literature . accounted for in the reconstruction algorithm. For the in vivo studies, the experimental data were acquired from a mouse trunk using a small animal imaging system that has been described in detail in earlier works . As opposed to the previous studies where the medium was assumed to be homogeneous, in this study the acoustic variation between the background and the bulk mouse tissue were accounted for in the reconstruction algorithm. A reconstructed initial pressure distribution image for fixed constant sound speed value is shown in Fig. 5(a) . In the reconstructed image, strong surface and interior vessel structures are observed. While the interior vessel structures appear out of focus, the surface vessels appear to be in focus. When assuming a homogeneous acoustic media, the image reconstruction algorithm could not produce images where both the surface and the interior vessels could be concurrently focused with a single tuned speed of sound value. To account for this, an image reconstruction algorithm that assigned different speed of sound values to the bulk tissue and the background was employed. The PLS-TV algorithm that was based on the imaging model defined in Eqn. (28) was employed to reconstruct the initial pressure distribution. The reconstructed initial pressure distribution when a heterogeneous speed of sound distribution is used is shown in Fig. 5 