Translator Disclaimer
31 July 2013 Least squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography
Author Affiliations +
A computationally efficient approach that computes the optimal regularization parameter for the Tikhonov-minimization scheme is developed for photoacoustic imaging. This approach is based on the least squares-QR decomposition which is a well-known dimensionality reduction technique for a large system of equations. It is shown that the proposed framework is effective in terms of quantitative and qualitative reconstructions of initial pressure distribution enabled via finding an optimal regularization parameter. The computational efficiency and performance of the proposed method are shown using a test case of numerical blood vessel phantom, where the initial pressure is exactly known for quantitative comparison.

Photoacoustic (PA) imaging is an emerging, noninvasive, in vivo biomedical imaging modality.1,2 A nanosecond laser pulse is generally used to irradiate biological tissue, resulting in a temperature rise from optical absorption and subsequently producing pressure waves due to thermoelastic expansion.2 The pressure waves are then acquired using a wide-band ultrasonic transducer at various locations around the surface of the tissue. A reconstruction algorithm is deployed that maps the initial pressure rise (proportional to the absorbed optical energy density) within the tissue from the recorded PA signals.3

Several PA image reconstruction algorithms were proposed in the literature,4,5 including analytical algorithms in the form of filtered back projection (BP) or algorithm based on Fourier transform. Their limitations include the requirement of large amount of data and are limited in terms of quantitative estimation.68 To overcome this limitation, various iterative image reconstruction algorithms have been proposed68 to improve the quantative accuracy of the reconstructed images, at the same time being computationally efficient. Moreover, in case of full-view data sets, the least squares QR (LSQR)-based reconstruction scheme was used, which indirectly provides regularized solution with an added advantage of being highly efficient.7 The limited-view data is reconstructed using a standard Tikhonov regularization, which is time consuming and requires an explicit regularization parameter.58

In this letter, we propose a Tikhonov regularization framework based on LSQR decomposition, where Q and R represent an orthogonal and upper triangular matrices, respectively, which uses Lanczos bidiagonalization to provide dimensionality reduction to the system of equations in the case of limited-view data set. This is further used to carry out a simplex method-based optimization procedure to find the optimal regularization parameter. The performance of the proposed method is compared with the generalized cross validation (GCV)9 and L-curve9 methods, along with the analytical methods, such as BP4 and k-wave-based time-reversal reconstruction10 using a numerical blood vessel network phantom.

The system matrix approach has been adopted here to describe the PA data collection process, which can be represented by a Toeplitz matrix of a time-varying causal system. The image (dimension of n×n) is converted into a long vector by stacking all columns one below the other, represented by x (dimension of n2×1). The system matrix (A) has a dimension of m×n2. Here, each column of A represents the impulse response corresponding to each pixel in the image. Moreover, the time-varying data is stacked to result in a long vector having dimensions m×1 (which makes the number of rows of A to be m). In order to improve the computation time for building the system matrix, the system response was measured only once for the corner pixel [making x(corner pixel)=1 and the rest of the entries were made zero], which forms the first column entry of A. The rest of the columns (n21) are filled by using shifting and attenuation properties of the PA signal. This approach assumes that the medium has homogeneous ultrasound properties.

The system response for the corner pixel is recorded using k-wave matlab toolbox,10 which simulates the PA wave prorogation in two dimension. The simulation geometry had a computational grid of 701×701pixels (0.1mm/pixel). Forty detectors were placed in a circular fashion of 34-mm radius. Each detector was assumed to be a point detector with a frequency response of 2.25 MHz as center frequency and 70% bandwidth. The imaging region was restricted to 201×201pixels located at the center, resulting in n2=40,401. A time step of 50 ns having 1000 time steps was used in recording each signal (making m=40,000). The simulations assumed a sound speed of 1500m/s.

In summary, the forward model of PA imaging can be written as

Eq. (1)

where x is a long column vector (unknown, representing the initial pressure rise, p0) and b is a measurement vector. The simple BP (analytical) image reconstruction scheme becomes4

Eq. (2)

where T represents the transpose of the matrix. As it is noniterative, this method is computationally efficient but known to provide only qualitative results.4

Both BP and time-reversal methods are analytical in nature but lack the quantitative nature of the results.5 In cases of limited data, typically a model-based reconstruction is employed which relies on minimizing the data-model misfit along with a regularization function, therefore the objective (cost) function in this case can be written as

Eq. (3)

where λ is the regularization parameter. The 2-norm is represented by .22. The function Ω is minimized with respect to x, leading to a direct solution

Eq. (4)


The GCV method9 is the most popular automated approach for estimating the optimal regularization parameter λopt using following function

Eq. (5)

where the SVD of A=HΣGT where, Σ is a diagonal matrix containing singular values (σ). The left and right orthogonal matrices are given by H and G, respectively. The L-curve is another popular scheme for estimating the optimal regularization parameter.9 The corner of L-curve gives the optimal regularization value, and as with GCV method, it does not require any prior information. The LSQR method is one of the variants of the conjugate gradient method used to solve a large system of equations. One of the main contributions of this letter is to use the LSQR-type algorithm to optimally determine the regularization parameter in PA image reconstruction. This is accomplished by using a Lanczos bidiagonalization of the system matrix (A). The left and right Lanczos matrices along with the bidiagonal matrix are related to the system matrix as shown below:11,12

Eq. (6)


Eq. (7)

where B represents the lower bidiagonal matrix, and U and V represent the left and right orthogonal Lanczos matrices, respectively. The unit vector of dimension k×1 is represented by ek (=1 at the k’th row and 0 elsewhere). Note that the dimensions of Uk and Vk are (m×k) and (n2×k), with k representing the number of iterations in the bidiagonalization procedure. Finally, ui and vi represent the left and right Lanczos vectors. The Bk is the bidiagonal matrix having α1,,αk in the main diagonal and β1,,βk in the lower subdiagonal of the matrix with a dimension of [(k+1)×k].

Now the Tikhonov minimization update for the equation for the LSQR-type method12 is given by [which is equivalent to Eq. (4)]:

Eq. (8)


Here β0 is the 2-norm of b. Once, xk (reduced x) is estimated then the initial pressure can be obtained using the relation p0=x=Vkxk.

Determination of the optimal number of Lanczos iterations (kopt) and optimal regularization parameter λopt is given in Algorithm 1. The advantage of LSQR-type method in finding the initial pressure rise distribution p0 lies in its dimensionality reduction capability which makes the update as xk [Eq. (8)] with kn2. The major role in the entire optimization procedure is characterized by k (number of Lanczos bidiagonalization). This factor determines the size of the bidiagonal matrix, Bk [dimension of (k+1)×k].11,12 In this work, the optimal number of Lanczos iterations turned out to be 25 obtained using Algorithm 1, making kopt=25. The optimal λ is searched within the specified bound (λlim) and is set to 1000. A gradient-free simplex method type algorithm is used due to its computational efficiency to compute the optimal regularization parameter (λopt).12 The λopt is found corresponding to k=25 is chosen for LSQR method. Using these in Eq. (8) gave xk and subsequently p0 (x).

Algorithm 1

Algorithm for determining optimal number of Lanczos iterations and optimal regularization.

• Input: Lanczos Bidiagonal Matrix Bk; Vk (k=1,2,,50); b; β0; A; λlim.
• Output: Optimal number of Lanczos iterations: kopt and optimal regularization parameter: λopt
• for k=1,2,,50 do Steps 1–3
1. Estimate the optimal λ for the given k(λkopt) Simplex method is used to find λkopt in the range of [0 λlim], with x=Vk*xk, found using Eq. (8).
2. Compute xk with λ=λkopt using Eq. (8). x=Vkxk
3. Estimate resk=bA*x22kopt=index of minimal value ofresk and λopt=λkoptopt

A Linux workstation with Intel Xeon Dual Quad Core 2.33 GHz processor having 64 GB memory was used in all computations carried out in this work.

In order to show the effectiveness of the proposed method, a numerical blood vessel phantom was chosen. Note that measuring actual p0 in the experimental phantom case is extremely challenging, which makes the comparison of performance of the methods discussed here difficult. Figure 1(a) shows the blood vessel network used as a numerical phantom with a maximum initial pressure rise of 1 kPa. The k-wave tool box10 was used to generate the simulated PA data with 40 detectors around the object of interest. Subsequently, the simulated data had a signal-to-noise-ratio of 40 dB (1% noise) to mimic the real experimental situation. The reconstruction results obtained using various methods are shown in Fig. 1(b)1(f). The proposed method’s result is shown in Fig. 1(f), with the value obtained for λopt shown in the parenthesis at the top of the image. Figure 1(g) shows the one-dimensional cross-sectional plot for the reconstructed PA image using GCV, L-curve, and LSQR-based methods as well as target image [along the dotted line in Fig. 1(a)], quantitatively showing an improvement of at least 10 times in the recovered p0. As seen from Fig. 1, the performance of LSQR-type method is superior in terms of quantitation compared to their counterparts. The total computational time recorded for all reconstruction methods shown in Fig. 1 are 129.8, 1.3, 7516.5, 7130.7, 444.9 s, respectively (system matrix building time = 181 s). A speed up factor of 17 was achieved by the proposed method compared to GCV method. The results indicate that among the model-based methods, the proposed method (LSQR) is the most efficient and promising technique for real-time imaging.

Fig. 1

(a) Numerical blood vessel phantom is used for the study (dimensions are in measured in millimeters). (b–f) Reconstructed photoacoustic images using k-wave interpolated, backprojection, reconstructions using λ [Eq. (4)] obtained by GCV [Eq. (5)], L-curve and LSQR [Eq. (8)] methods, respectively. (g) One-dimensional cross-sectional plot for the results presented in (a),(d),(e), and (f) along the dotted line shown in (a).


It is important to note that the LSQR (without explicit regularization) has been extensively used in full-view data cases (where mn2) and shown to be effective in quantitation.68 In limited data cases (where mn2), the Tikhonov minimization scheme is more effective with the important condition of finding an optimal regularization parameter. This work addresses this problem of finding an optimal regularization parameter automatically without prior information with an added advantage of reducing the dimensionality of the problem, making it a highly computationally efficient method.



L. H. V. WangS. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science, 335 (6075), 1458 –1462 (2012). SCIEAS 0036-8075 Google Scholar


M. Pramaniket al., “Design and evaluation of a novel breast cancer detection system combining both thermoacoustic (TA) and photoacoustic (PA) tomography,” Med. Phys., 35 (6), 2218 –2223 (2008). MPHYA6 0094-2405 Google Scholar


J. Gamelinet al., “Curved array photoacoustic tomographic system for small animal imaging,” J. Biomed. Opt., 13 (2), 024007 (2008). JBOPFO 1083-3668 Google Scholar


G. Paltaufet al., “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acoust. Soc. Am., 112 (4), 1536 –1544 (2002). JASMAN 0001-4966 Google Scholar


K. Wanget al., “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Biol., 57 (17), 5399 –5424 (2012). PHMBA7 0031-9155 Google Scholar


X. L. Dean-BenV. NtziachristosD. Razansky, “Acceleration of optoacoustic model-based reconstruction using angular image discretization,” IEEE Trans. Med. Imag., 31 (5), 1154 –1162 (2012). ITMID4 0278-0062 Google Scholar


X. L. Dean-Benet al., “Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imag., 31 (10), 1922 –1928 (2012). ITMID4 0278-0062 Google Scholar


A. Buehleret al., “Model-based optoacoustic inversions with incomplete projection data,” Med. Phys., 38 (3), 1694 –1704 (2011). MPHYA6 0094-2405 Google Scholar


R. P. K. JagannathP. K. Yalavarthy, “Minimal residual method provides optimal regularization parameter for diffuse optical tomography,” J. Biomed. Opt., 17 (10), 106015 (2012). JBOPFO 1083-3668 Google Scholar


B. E. TreebyB. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt., 15 (2), 021314 (2010). JBOPFO 1083-3668 Google Scholar


C. C. PaigeM. A. Saunders, “LSQR: An algorithm for sparse linear equations and sparse least squares,” ACM Trans. Math. Softw., 8 (1), 43 –71 (1982). ACMSCU 0098-3500 Google Scholar


J. PrakashP. K. Yalavarthy, “A LSQR-type method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography,” Med. Phys., 40 (3), 033101 (2013). MPHYA6 0094-2405 Google Scholar
© 2013 Society of Photo-Optical Instrumentation Engineers (SPIE) 0091-3286/2013/$25.00 © 2013 SPIE
Calvin B. Shaw, Jaya Prakash, Manojit Pramanik, and Phaneendra K. Yalavarthy "Least squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography," Journal of Biomedical Optics 18(8), 080501 (31 July 2013).
Published: 31 July 2013

Back to Top