Accelerating cross-validation with total variation and its application to super-resolution imaging

We develop an approximation formula for the cross-validation error (CVE) of a sparse linear regression penalized by ℓ1-norm and total variation terms, which is based on a perturbative expansion utilizing the largeness of both the data dimensionality and the model. The developed formula allows us to reduce the necessary computational cost of the CVE evaluation significantly. The practicality of the formula is tested through application to simulated black-hole image reconstruction on the event-horizon scale with super resolution. The results demonstrate that our approximation reproduces the CVE values obtained via literally conducted cross-validation with reasonably good precision.


Introduction
At present, in many practical situations of science and technology, large high-dimensional observational datasets are created and accumulated on a continuous basis. An essential difficulty concerning the treatment of such high-dimensional data is the extraction of meaningful information. Sparse modeling [1,2] is a promising framework for overcoming this difficulty, which has recently been utilized in many disciplines [3,4]. In this framework, a statistical or machine-learning model with a large number of parameters (explanatory variables) is fitted to the data, in conjunction with a certain sparsity-inducing penalty. This penalty should be appropriately chosen with consideration of the processed data. One representative penalty is the ℓ 1 regularization, which retains certain preferred properties, such as the statistical model convexity [5,6]. A similar penalty that has received more recent focus is the so-called "total variation (TV)" [7][8][9], which can be regarded as the ℓ 1 regularization imposed on the difference between neighboring explanatory variables. The TV yields "continuity" of the neighboring variables, which is suitable for the processing of certain datasets expected to have such continuity, such as natural images [4,[7][8][9].
Another common difficulty associated with the use of statistical models is model selection. In the context of image processing using the ℓ 1  during the selection of appropriate regularization weights. A practical framework to select these weights, which is applicable to general situations, is cross-validation (CV). CV provides an estimator of the statistical-model generalization error, i.e., the CV error (CVE), using the data under control, and the minimum CVE obtained when sweeping the weights yields the optimal weight values. This versatile framework is, however, computationally demanding for large datasets/models, and this problem frequently becomes a bottleneck affecting model selection. Thus, reducing the CVE computational cost could have a significant impact on a broad range of sparse modeling applications in various disciplines. Considering these circumstances, in this paper, we provide a CVE approximation formula for a statistical model of linear regression penalized by the ℓ 1 and TV terms, to efficiently reduce the computational cost. Note that the formula for the case penalized by the ℓ 1 term alone has already been proposed in [10], and the formula presented herein is a generalization of it. Below, we show the formula derivation and perform a demonstration in the context of super-resolution imaging. The processed images employed in this study are reconstructed from simulated observations of black holes on the event-horizon scale for the Event Horizon Telescope (EHT, see [11][12][13]) full array. Note that our formula will be applied to actual EHT observations to be conducted after April 2017.

Problem setting
Let us suppose that our measurement is a linear process, and denote the measurement result as y 2 R M and the measurement matrix as The explanatory variables, corresponding to the images that will be examined in the later demonstration, are denoted by x 2 R N . The quality of the fit to the data is described by the residual sum of squares (RSS), i.e., Eðxjy; AÞ ¼ 1 2 jjy À Axjj 2 2 . In addition, we consider the following penalty consisting of ℓ 1 and TV terms: where the T(x) term corresponds to the TV and is expressed as and @i denotes the neighboring variables of the ith variable. There is some variation in the definition of "neighbors"; here, we follow the standard approach [7][8][9]. That is, x is assumed to be a two-dimensional image and the neighbors of the ith pixel correspond to the right and down pixels. However, the bottom row (the rightmost column) of the image is exceptional, as the neighbor of each pixel in that row (column) corresponds to the right (down) pixel only. Note that the developed approximation formula presented below is independent of this specific choice of neighbors and can be applied to general cases. For this setup, we consider the following linear regression problem with the penalty given in Eq (1)x where arg min u ff ðuÞg generally represents the argument that minimizes an arbitrary function f(u). Further, we consider the leave-one-out (LOO) CV of Eq (3) in the form Note that the system without the μth row of y and A is referred to as the "μth LOO system" hereafter. In this procedure, the CVE, i.e., the generalization error estimator, is where a > m ¼ ðA m1 ; Á Á Á ; A mN Þ is the μth row vector of A. We term this simply the "LOO error (LOOE)." Computing the LOOE requires solution of Eq (4) M times, by definition, which is computationally expensive. Therefore, the purpose of this paper is to avoid this computational expense by deriving an approximation formula of Eq (5).

Approximation formula for softened system
When M is sufficiently large, i.e., the number of observations is large enourgh, the difference between the LOO solutionx nm and the full solutionx is expected to be small. This intuition naturally motivates us to conduct a perturbation connecting these two solutions. To conduct this perturbation, we "soften" the penalty by introducing a small cutoff δ(> 0) in the TV, having the form where An approximation formula in the presence of ℓ 1 regularization with smooth cost functions has already been proposed in [10]. We employ that formula here and take the limit δ ! 0.
To state the approximation formula, we begin by defining "active" and "killed" variables. Owing to the ℓ 1 term, some variables are set to zero inx; we refer to these variables as "killed variables." The remaining finite variables are termed "active variables." We denote the index sets of the active and killed variables by S A and S K , respectively. The active (killed) components of a vector x are formally expressed as x S A (x S K ). For any matrix X, we use double subscripts in the same manner. For example, for an N × N matrix, a submatrix having row and column components of S A and S K , respectively, is denoted by X S A SK .
The approximation formula can be derived through the following two steps. Note that, in this derivation, a crucial assumption is that the sets of active and killed variables are common among the full and LOO systems. This assumption may not hold exactly in practice, but the resultant formula is asymptotically exact in the large-N limit [10].
The first step is to compute the values of the active variables and their response to small perturbation. The active variables are determined by the extremization condition of the softened cost function with respect to the active variables, such that The focus here is the response of this solution when a small perturbation −h Á x is incorporated into the cost function. A simple computation demonstrates that the active-active components of the response function, , are equivalent to the inverse of the costfunction Hessian where @ 2 x denotes the Hessian operator @ 2 x @ 2 @x i @x j and G nμ is the Gram matrix of A nμ , i.e., The other components of the response function are identically zero, from the stability assumption of S K and because the killed variables are zero, withx S K ¼x nm S K ¼ 0. In the second step, we connect the full solution to the LOO solution, through the above perturbation with an appropriate h. To specify the perturbation, we assume that the difference d d ¼x d Àx dnm is small and expand the RSS of the full system with respect to d δ as follows: This equation implies that the perturbation between the full and LOO systems can be expressed as h m ¼ ðy m À a > mx d Þa m . Hence, we obtain The Hessian of the full system has a simple relationship with the LOO Hessian, such that where the approximation at the righthand side comes from replacingx d withx dnm in the argument of R δ (x). Inserting Eqs (12 and 13) in conjunction with w d using the Sherman-Morrison formula for matrix inversion, we find According to Eq (14), we can compute the LOOE only from the full solutionx d , without actually performing CV, which facilitates considerable reduction of the computational cost.

Handling a singularity
Let us generalize Eq (14) to the limit δ ! 0, where the penalty contains another singular term in addition to the ℓ 1 term. This TV singularity tends to "lock" some of the neighboring variables, i.e., x j = x i (8j 2 @i), which corresponds to t i = 0 in Eq (2). If two different vanishing TV terms, t i and t j , share a common variable x r , all the variables in those TV terms take the same value ). In this manner, the active variables are separated into several "locked" clusters, with all the variables inside a cluster having an identical value. This implies that the variable response to a perturbation, χ = lim δ!0 χ δ , should have the same value for all variables in a cluster and may, therefore, be merged. Below, we demonstrate this behavior for the δ ! 0 limit. For the derivation, we assume that the clusters are common to both the full and LOO systems, similar to the assumption for S A and S K . For convenience, we index the clusters by α, β 2 C and denote the number of clusters by |C|; the index set of variables in a cluster α is represented by S α and the total set of indices in all clusters is denoted by S C [ α S α . Hereafter, we concentrate on the active variable space only and omit the killed variable space. The complement set of S C , i.e., the set of isolated variables that do not belong to any cluster, is denoted by S I and, thus, Two crucial observations for the derivation are the "scale separation" and the presence of the "zero mode." For vanishing TV terms, a natural scaling to satisfy lim Once this scaling is assumed, we realize that the components of the Hessian that are directly related to the clusters diverge. Let us define byŜ a the set of TV terms corresponding to cluster α, i.e.,Ŝ a ¼ fijðfig [ @iÞ & S a g. Hence, by construction and for all and F δ consists of the remaining O(1) terms. This decomposition can be schematically expressed as We denote the basis of the current expression by {e i } i2S A , with (e i ) j = δ ij , and move to another basis that diagonalizes D d S C S C . Each D d a has a "zero mode," and its normalized eigenvector is given by for i 2 S α and 0 otherwise, in the full space. This behavior originates from the symmetry, such that the ft d i g i2Ŝ a are invariant under a uniform shift in the S α sub-space, i.e., x j ! x j + Δ (8j 2 S α ) for 8Δ 2 R. This invariance can also be directly seen from a property of the Hessian, i.e., @ 2 In addition, we represent the set of normalized eigenvectors of all the other modes of D d a , which have eigenvalues λ αa that are proportional to 1/δ and positively divergent, as fu aa g jS a jÀ 1 where relevant for the evaluation of (H δ ) −1 . These considerations yield the explicit formula of χ as As the non-zero components of the zero mode z α are identically given as 1= ffiffiffiffiffiffiffi jS a j p , all these coefficients can be easily expressed by the original coefficients F ij , as and F iα = F αi by the symmetry. Now, all the components are explicitly specified. The form of χ in the original basis {e i } i2S A can be accordingly assessed by moving back from the basis {{{u αa } a , z α } α , {e i } i2S I }, which completes the computation. Some additional consideration of the above computation demonstrates that we can shorten some steps and obtain a more interpretable result. We introduce a jS IþZ j Â jS IþZ j matrix " F as with the remaining components being identical to those of FS IþZSIþZ , i.e., " F S I S I ¼ F S I S I . Eqs (19) and (20) indicate that " F is simply the matrix summing the rows and columns in each cluster to a row and a column. It is natural that " F has a direct connection to χ, because the locked variables in a cluster should exhibit the same response against perturbation. In fact, the response function χ in the original basis is expressed using " F as This can be directly shown from Eqs (17 and 19), using the relation FS g a gÞ, and the blockwise matrix inversion formula. Eqs (14) and (21) constitute the main result of this paper.

Numerical stability and the softening constant δ
For handling the singularity of the cost-function Hessian, we have introduced the softening constant δ in the TV and finally taken the δ ! 0 limit. In practical implementations, however, we should keep δ small but finite. To see the reason, it is sufficient to see a simple example with just three variables fx i g 3 i¼1 . The softened TV is defined as where p = x 2 − x 1 , q = x 3 − x 1 are introduced. The corresponding gradient and Hessian are The zero point of the gradient is given by p = q = 0 irrespectively of the δ value. Inserting this into the Hessian, we get one zero mode proportional to (1, 1, 1) > and two finite modes whose eigenvalues are (3/δ, 1/δ) being divergent in the δ ! 0 limit. This exactly matches with the assumptions of the approximation formula.
On the other hand, if we first take the limit δ ! 0 before taking the zero gradient limit p, q ! 0, we see that two zero modes appear: One is proportional to (1, 1, 1) > and the other is to (p + q, q − 2p, p − 2q) > . This is a bad news because the second zero mode, which remains even in the limit p, q ! 0, is never taken into account when deriving the approximation formula: The derivation essentially depends on how the zero mode behaves and our formula loses its justification if such unexpected zero modes exist.
These considerations manifest that the two limits, lim δ!0 and lim p,q!0 , are not exchangeable in the TV Hessian. The derivation of our approximation formula assumes lim δ!0 lim p,q!0 and thus the algorithmic implementation should reflect this limit in a certain way. A simple way is to keep δ small but finite, which is actually a common technique to enhance the numerical stability when using the TV [14]. The choice of the amplitude of δ is related to the numerical precision when solving the optimization problem (3). A practical choice is stated in the next subsection.

Procedures
Here, we state the procedures for implementation of Eqs (14 and 21) in a numerical computation. Suppose that we have an algorithm to solve Eq (3) and to provide the solutionx given y, A, λ ℓ 1 , and λ T . Using this solution and introducing a finite δ in the Hessian by the reason discussed above, we can assess the LOOE through the following steps: 1. The sets of active and killed variables, S A and S K , are specified fromx.

A new index set S R = {{α} α2C
, S I } is defined.
6. On S R , the merged Hessian " F is constructed from F, as " Similarly, the merged measurement matrix " A is defined as " 7. Using " F and " A, the LOOE factor in Eq (14) is computed as where " a T m is the μth row vector of " A and x = A\b is the solution of the linear equation Ax = b. 8. Using the LOOE factor andx, the LOOE is evaluated from Eq (14).
At step 7, we take the left division " F S R S R n" a mS R instead of the inverse w ¼ " F À 1 for numerical stability. The cluster enumeration at step 3 involves a delicate point in the definition of C and {S α } α2C . Because of the limited precision in the numerics, the TV term jx j Àx i j ðj 2 @iÞ never exactly vanishes; therefore, we need a certain threshold to extract the cluster structure from the TV terms. Here, we introduce the threshold θ and enumerate the clusters as follows: 3-2. An empty set C = ϕ is prepared and the cluster index α = 1 is defined.
3-3. The following steps are repeatedly implemented while L is non-empty: (i). Two empty sets, S tmp = ϕ and S cluster = ϕ, are prepared; (ii). One link is selected and removed from L. The variable indices in the link are entered into S tmp ; (iii). The following steps are repeatedly implemented while S tmp is non-empty: a. One index i in S tmp is selected and moved from S tmp to S cluster ; b. If the above chosen index i exists in S L , all the links to i are removed from L, and S L is updated accordingly. The variables linked to i are entered into S tmp ; c. S tmp S tmp − S cluster .
(iv). The variables in S cluster constitute a cluster. S α = S cluster is defined and α is entered into C; (v). α α + 1.
The entire procedure presented above implements Eqs (14 and 21). A debatable point would be the values of θ and δ. In most of iterative algorithms as the one in [8,9], there is an inevitable finite error of the TV term even when it should vanish. Let us express the "scale" of this error as t i ðxÞ % y num > 0. By construction, the threshold θ is related to this numerical error and it is appropriate to choose θ % θ num ; the softening constant δ should be sufficiently larger than θ num because it does implement the assumed order of two limits, lim δ!0 lim p,q!0 , in derivation of the approximation formula. Overall, the relation must be satisfied. We have numerically checked how strict this principled relation is, and found that the approximation result is not sensitive to the choice of θ as long as it is sufficiently smaller than δ. Although a little more delicate points are involved in the choice of δ, we have also found that in a wide range of δ the approximation result is stable and the cost-function Hessian is safely invertible. Based on these observations, in the application of our formula below, the default values are set to be δ = 10 −4 and θ = 10 −12 . They are chosen according to our datasets and experimental setup: The maximum value of the non-softened TV terms is scaled as max i t i ðxÞ ≳ 10 À 4 and the numerical precision is about θ num % 10 −12 ; the former value is reflected to δ and the latter one is used in θ. Coincidently, this default value of δ accords with the one in [14]. The examination result of the sensitivity to δ and θ will be reported below.
Another noteworthy point is that these procedures can be easily extended to other variants of the TV. For example, for the so-called anisotropic TV [9], T ani = ∑ i ∑ j2@i |x j − x i |, we set F = G in step 4 and modify the definition of the link in step 3-1 accordingly, so as to render our formula applicable. In the case of the square TV, T sq = ∑ i ∑ j2@i (x j − x i ) 2 (1/2)x > J x, the formula can be significantly simpler, because this TV has no sparsifying effect and the formula of the simple ℓ 1 case can be employed. We can employ Eq (14) with χ S A S A = (G S A S A + λ T J S A S A ) −1 directly, without the need for cluster enumeration.

Application to super-resolution imaging
To test the usefulness of the developed formula, let us apply the derived expression to the super-resolution reconstruction of astronomical images. A number of recent studies have demonstrated that sparse modeling is an effective means of reconstructing astronomical images obtained through radio interferometric observations [15][16][17][18]. In particular, the capability of high-fidelity imaging in super-resolution regimes has been shown, which renders this technique a useful choice for the imaging of black holes with the EHT [17][18][19][20][21]. We adopt the same problem setting as [17,20] and demonstrate the efficacy of our approximation formula through comparison with the literally conducted 10-fold CV result. Here, x i denotes the ith pixel value and A is (part of) the Fourier matrix. The dataset y is generated through the linear process where ξ is a noise vector and x 0 is the simulated image, which we infer given y and A.
In this work, we use data for simulated EHT observations based on three different astronomical images, which are available as sample data for the EHT Imaging Challenge. Our datasets 1, 2, and 3 correspond to the sample datasets 1, 2, and 5, respectively, available from [22] at July 2017. The images are reconstructed with N = 10000 = 100 × 100 pixels and with 160, 250, and 100 μ as fields of view, which are identical to the original images of Datasets 1, 2, and 3 from the EHT Imaging Challenge, respectively. We test four values for each λ ℓ 1 and λ T : λ ℓ 1 2 (M/2) × {1, 10, 100, 1000} and λ T 2 (M/8) × {1, 10, 100, 1000}. M is 1910, 1910, and 2786, for Datasets 1-3, respectively. Later, we also use different size data from the same datasets, for checking the size dependence of the result. Table 1 shows the mean CVE values for the three datasets, determined by the 10-fold CV and by our approximation formula for varying λ T . λ ℓ 1 is fixed to the optimal value, which is coincidently common for all datasets and satisfies 2λ ℓ 1 /M = 1. It is clear that the approximate CVE values accord well with the 10-fold results, even on the error-bar scale, demonstrating that our approximation formula works very well. Note that the error bar for the approximation is given by the standard deviation of the M terms in Eq (14) divided by ffiffiffiffiffiffiffiffiffiffiffiffiffi M À 1 p . To directly observe the reconstruction quality, in Fig 1 we display the images at all investigated parameters and the reconstructed image at the optimal λ ℓ 1 and λ T for Dataset 3, as well as the associated errors plotted against λ ℓ 1 and λ T in Fig 2. Again, we can see the proposed method approximates the 10-fold result well, and the reconstructed image reasonably resembles the original. The RSS is monotonic with respect to the changes of λ l 1 and λ T but the approximate LOOE is not, which implies that the LOOE factor computed through Eq (21) appropriately reflects the effect of the penalty terms.
Next, we check the sensitivity of the approximate result to the tuning constants δ and θ. In Fig 3, the approximate LOOEs at the optimal λ ℓ 1 are plotted against λ T when changing δ (left) and θ (right). This indicates that the approximate LOOEs are stable against the change of both δ and θ. Hence, we may choose these values rather arbitrarily. This is a good news because tuning them makes the problem more numerically amenable: Enlarging δ makes the computation Table 1. CVE values determined by 10-fold CV and our approximation formula against λ T . λ ℓ 1 is fixed to the optimal value (2λ ℓ 1 /M = 1, coincidentally common to all cases). The number in brackets denotes the error bar to the last digits. The optimal values are bolded. The tuning constants δ and θ are set to be δ = 10 −4 and θ = 10 −12 , respectively.   1, 10)). The images are convolved with a circular Gaussian beam on the right-hand side, the full width at half maximum (FWHM) of which is 25% of the nominal angular resolution of the EHT and corresponds to the diameters of the yellow circles. This coincides with the optimal resolution minimizing the mean square error between them.
of the Hessian inversion more numerically stable; increasing θ lowers the effective degrees of freedom. The second property associated with θ is really beneficial when treating a large-size dataset, because it can downsize the Hessian and reduce the cost for computing its matrix inversion. In Table 2, the values of the effective degrees of freedom are given when changing θ.
The reduction of the degree of freedom at large (yet small enough compared to δ = 10 −4 ) θ is significant, which encourages us to apply the proposed formula to larger-size datasets. Finally, let us see the data-size dependence of the approximation accuracy and of the computational cost for solving Eq (3) and for obtaining the approximate LOOE from the solution. The data analyzed here is an identical simulated image of black hole expressed with different number of pixels. When solving Eq (3), we used Intel(R) Core(TM) i7-5820K CPU of 3.30GHz with 6 cores for N = 50 2 = 2500 and Intel(R) Xeon(R) CPU E5-2699 v3 of 2.30GHz with 36 cores for N = 100 2 and 150 2 , and employed an algorithm called "MFISTA" proposed in [8,9]. Meanwhile, we used a laptop of a 1.7 GHz Intel Core i7 with two CPUs for evaluating  the approximate LOOE. Hence the comparison is not fair and unfavorable to the approximation formula. The left panel indicates that the approximation accuracy becomes better for larger sizes. This is reasonable because the perturbation we have employed should have better accuracy as the model and data become larger, though the accuracy at N = 50 2 = 2500 is already good. The right panel clearly shows the advantage of the developed formula: The actual computational time of the approximate LOOE is significantly shorter than that of the algorithm convergence for solving Eq (3) in the investigated range of system sizes, even under the unfair comparison mentioned above. However, this advantage will be less prominent if the model becomes very large: Our approximation formula needs the Hessian inversion whose computational cost is scaled as O((|C| + |S I |) 3 ) % O(N 3 ), while MFISTA requires the cost of O(N 2 ) as long as the number of steps to convergence is constant against N. The crossover size at which these two computational costs become comparable is roughly estimated as N × % 10 6 , though such crossover tendency cannot be seen yet from Fig 4. For such large systems, a new fundamental solution should be tailored to resolve the computational-cost problem, though tuning θ to a large value in the present method can still be a good first aide.

Conclusion
In this paper, we have developed an approximation formula for the CVE of a sparse linear regression penalized by ℓ 1 and TV terms, and demonstrated its usefulness in the reconstruction of simulated black hole images. Our derivation is based on the perturbation assuming the small difference between the full and leave-one-out solutions. This assumption will not be fulfilled for some specific cases, i.e. when the measurement matrix is sparse. However, for most of dense measurement matrices, such as the Fourier matrix discussed in this paper, our Table 2. The effective degrees of freedom jS IþZ j, the number of clusters + the number of isolated variables, against θ for Dataset 3 at δ = 10 −4 and the optimal parameters (2λ ℓ 1 , 8λ T )/M = (1, 10).  assumption will be reasonably satisfied. Hence we expect the range of application of our formula is wide enough and we would like encourage the readers to use this formula in their own work. It is also straightforward to generalize the developed formula to other types of TV, and two examples of the generalization for the anisotropic and square TVs have been explained.
The key concept of our formula, perturbation between the LOO and full systems, is very general and can be applied to more general statistical models and inference frameworks [23]. The development of practical formulas for those cases will facilitate higher levels of modeling and computation.