{'fraction_non_alphanumeric': 0.04811890872253467, 'fraction_numerical': 0.015474901316938555, 'mean_word_length': 4.0210470928702975, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning multiple regularization parameters for generalized Tikhonov regularization using multiple data sets without true data

During the inversion of discrete linear systems, noise in data can be amplified and result in meaningless solutions. To combat this effect, characteristics of solutions that are considered desirable are mathematically implemented during inversion, which is a process called regularization. The influence of provided prior information is controlled by non-negative regularization parameter(s). There are a number of methods used to select appropriate regularization parameters, as well as a number of methods used for inversion. In this paper, we consider the unbiased risk estimator, generalized cross validation, and the discrepancy principle as the means of selecting regularization parameters. When multiple data sets describing the same physical phenomena are available, the use of multiple regularization parameters can enhance results. Here we demonstrate that it is possible to learn multiple parameter regularization parameters using regularization parameter estimators that are modified to handle multiple parameters and multiple data. The results demonstrate that these modified methods, which do not require the use of true data for learning regularization parameters, are effective and efficient, and perform comparably to methods based on true data for learning the relevant parameters.


Introduction
We consider solutions of linear problems described by where A ∈ R m×n with m ≥ n and d is known, but represents d = b + η, with η being a realization of a random vector and b = Ax true . Even for invertible square matrices A, direct matrix inversion when A is ill-conditioned is not recommended due to the noise in the data. Rather, an alternative approach is regularization, in which desired characteristics of a solution are described mathematically and incorporated into the formulation with the aim to produce a more well-posed problem. The Tikhonov regularized solution, x(α), is defined as the solution of the Tikhonov regularization problem, [48], as given by Here, the scalar α > 0 is a regularization parameter and L is a q × n matrix representation of a linear operator. The term Lx 2 2 is an example of a penalty function [50], and L is called the penalty matrix. If L = I n , the n × n identity matrix, then regularization via eq. (2) is called standard or zeroth-order Tikhonov regularization [3]. Other standard choices of L include approximations of first and second order derivative operators [40,46,50]. Typically, q < n for one-dimensional (1D) problems without specific attention to boundary conditions, though our formulations throughout this paper do not assume we are only dealing with 1D problems.
The quality of x(α) depends on the choice of both α and L. There are a number of methods for selecting α when L has been fixed. As an example, the Morozov discrepancy principle (MDP) [38], which requires knowledge of the variance of the noise distribution of the data, selects the regularization parameter as the root of a function. The unbiased predictive risk estimator (UPRE), [33], which also requires knowledge of the variance of the noise distribution of the data, yields the regularization parameter as the minimizer of a function. The method of generalized cross validation (GCV), [51,52], which does not require the noise distribution be known, also yields the regularization parameter as a minimizer of a function. Other methods do not solve minimization or root-finding problems; for example, the L-curve method selects a regularization parameter as the value that locates the point of maximum curvature of a function [22,25].
There has also been a considerable amount of research in selecting sets of regularization parameters for a pre-selected set of penalty matrices, a process called multi-parameter (MP) regularization [7,11,19,32,54]. Approaches to MP regularization using versions of the L-curve and MDP methods can be found in [4] and [53], respectively. A MP GCV method was also considered in [36,37]. Windowing, either in the data domain or the frequency domain, can also be applied to determine multiple regularization parameters, windowing wavelet coefficients was considered in [16,44]. Examples of windowed regularization in other frequency domains, such as those generated by discrete trigonometric transforms or the singular value decomposition (SVD), have been presented in [10,13,28]. There is also recent work on learning seminorms as regularization operators [26].
Techniques for utilizing, and analyzing multiple data (MD) sets, permeate a multitude of scientific fields as diverse as geoscience [6,55], and the detection of cancers [43]. As noted, many of these techniques, including the UPRE, MDP, and GCV methods, have statistical foundations [27]. A comprehensive overview of data-driven approaches to inverse problems can be found in [2], while specific examples of applying multiple data sets for the solution of inverse problems include [1,9,11,21,29,47,49].
A main contribution of this paper is demonstrating and validating how the functions associated with the UPRE, MDP, and GCV methods can be adapted to handle multiple data sets, for both single (scalar) and multi-parameter regularization. These are chosen as representative methods that either require (UPRE and MDP), or do not require (GCV) prior knowledge of the statistics of the noise in the data. There is a distinction between the adapted methods and those that would use averages of the parameter selection functions. The development of the adapted methods is presented in section 3 and numerical results are shown in section 4. Sections 3.2 to 3.4 pertain to the UPRE, MDP, and GCV methods, respectively. In addition to developing these adapted parameter selection methods, results are presented regarding their relationship(s) with the original method(s) (i.e. the non-adapted methods), at the end of each subsection in section 3. Comparisons are also made to other regularization parameter selection methods that use multiple data sets; see section 4.
Most significantly, it is demonstrated through this work how multiple data sets can be used in conjunction with multi-parameter regularization using the windowing formulation introduced in [10]. Multi-parameter versions of the GCV method are presented in [10,37]; we extend these formulations for equivalent versions of the UPRE and MDP methods. Numerical results for the restoration of 1D and 2D signals, demonstrate that these new multiparameter regularization parameter estimators can be used for multiple data sets and that their performance competes with a learning approach in which the training stage requires the knowledge of true data, which is not the case here. Furthermore, parameters that have been obtained from one set of training images can be used on a separate set of validation (testing) images distinct from the original set, provided that the signal to noise ratios are close. Directions for future work are considered in section 5.

Summary of notation
Due to the consideration of multiple data sets as well as multiple parameter regularization for these data sets, the notation in this paper is extensive and will therefore be briefly summarized. We will use the acronyms "MP" and "MD" in the textual body of the paper and plots to refer to the multi-parameter and multi-data generalizations, respectively, of these methods and their corresponding functions. For example, the MD UPRE method means the extension of the UPRE method for multiple data sets. The multi-data multiparameter applications will be regarded using the combined acronym "MDMP". In the formulation, R denotes the number of data sets being considered, which are indexed by r. The index r may either be placed as a subscript or a superscript with parenthesis signs; for example, A (r) , x (r) , and d (r) represent the system matrix, solution, and data, respectively, associated with the rth system. The tilde (∼) is placed above matrices, vectors, or functions to denote their use for MD. Specifically, for matrices the tilde indicates the formation of a block diagonal matrix, e.g. A = diag(A (1) , . . . , A (R) ) for {A (r) } R j=1 . For vectors, the tilde indicates vertical concatenation, e.g. d is the ordered vertical concatenation of the vectors {d (r) } R j=1 . The functions used in each parameter selection method are indicated using F with subscripts denoting which method is being considered. For example, F UPRE (α) represents the scalar UPRE function for MD sets. P r represents the number of windows/parameters being considered for the rth data set; the individual windows are indexed by p. A general permutation matrix is represented by P . For matrices, vectors, and functions being used for MP regularization, the subscript "win" is added. Letters are bolded to indicate vectors; the term F UPRE win (α) represents the UPRE function for use in the ultimate case, where MD are being used for MP regularization with the parameters contained in the vector α.

Background
Here we first review the use of standard matrix decompositions for expressing the solution of the regularized problem eq. (2) for a single regularization parameter α. We then show how these same decompositions, combined with windowing techniques, can be used to provide compact expressions for the MP, and MDMP regularized solutions.

Singular Value Decomposition
The singular value decomposition (SVD) is a useful tool for analyzing the properties of methods for solving eq. (1), even if it is not used in practice for large scale problems, [20]. We consider the SVD of A in the full form given by where the m × m matrix U and the n × n matrix V are orthogonal. The columns of U and V denoted by U ·,j and V ·,j , respectively, are known as the left and right singular vectors of A. S is a m × n matrix with non zero entries s j on the leading diagonal ordered by s 1 ≥ s 2 ≥ . . . ≥ s r A ≥ 0. The s j are the singular values of A, and r A := rank(A) ≤ min{m, n}. If A is square (m = n) and invertible, then A −1 = V S −1 U T and the solution A −1 d can be written as where d = U T d is the vector of spectral coefficients of d.
When A is not square, or is square and not invertible, the pseudoinverse [41] A † = V S † U T can be used where S † ∈ R n×m is formed by reciprocating the nonzero diagonal elements of S T . When working with A † , it can be convenient to use the compact, or reduced, [3,30], form of the SVD given by where U r A ∈ R m×r A and V r A ∈ R n×r A consist of the first r A columns of U and V , respectively. Here S r A = diag(s 1 , . . . , s r A ) ∈ R r A ×r A is invertible. The columns of U r A form an orthonormal basis for range(A), while the columns of V r A form an orthonormal basis for range(A T ) [30, p. 340].
, the fundamental theorem of linear algebra [45] implies that the columns U 0 and V 0 form orthonormal bases for null(A T ) and null(A), respectively. Using eq. (4) the pseudoinverse of A is given by The summation representation of x † is the same as eq. (3) except that r A is the upper limit of summation.
In the case where null(A) = {0}, infinitely many solutions to Ax = d can be constructed as Generally, the summands in eq. (3) are numerically unstable for small s j when the coefficients d j do not decay as fast as the s j ; this describes the discrete Picard condition [23].
For an ill-conditioned matrix A, forming eq. (3) will often result in a meaningless solution. A common approach to overcome numerical instabilities is to multiply the summands in eq. (3) by filter functions φ j (α) that depend upon s j and a non-negative regularization parameter α. This provides the approximate solution where the matrix Φ(α) is diagonal with diagonal entries The most desired property of the filter functions is that φ j (α)/s j ≈ 1 for large values of s j and φ j (α)/s j ≈ 0 for small values of s j . A specific example is which is known as the standard Tikhonov filter function and corresponds to applying Tikhonov regularization to solve eq. (1), [48]. Equation (5), using filter functions given by eq. (6), is equivalent to the solution of the damped least squares problem [3]. Furthermore, eq. (7) solves the ordinary least squares problem where the system matrix has full rank for all α > 0. Using the SVD of A, the normal equations (A T A + α 2 I n )x = A T d corresponding to eq. (8) simplify to and eq. (5) is obtained by noting that (S T S + α 2 I n ) is a diagonal matrix that is guaranteed to be non-singular for α > 0. If A has full column rank and m ≥ n, then S T S + α 2 I n is nonsingular even when α = 0 but may be poorly conditioned. The invertibility of (S T S + α 2 I n ) provides a specific version of eq. (5), which is Analogous to eq. (8), the solution of eq. (2) for general L can be expressed in block form as If the matrix L in eq. (2) is non-singular, then the substitutions y = Lx, B = AL −1 , and c = d give y(α) = arg min This is known as the standard form of the regularization problem. Once y is obtained from eq. (11), the final solution is recovered by x = L −1 y. There are, however, many examples of matrices L that are singular, such as finite difference matrices that approximate derivative operators. In such cases, eq. (2) can still be transformed into the standard form eq. (11); this transformation was originally given by Eldén [17] and is discussed in [23].

Generalized Singular Value Decomposition
In light of having to consider two (often very different) matrices A and L in formulations such as eq. (2) and eq. (10), a brief discussion of the generalized singular valued decomposition (GSVD) is relevant. The discussion closely follows the presentation of the GSVD in [3], which uses the invertibility condition null(A) ∩ null(L) = {0}. This condition means that the block system matrix in eq. (10) has full column rank and x(α) is unique. Given real matrices A and L of size m × n and q × n, respectively, the decompositions exist where U is an m × m orthogonal matrix, V is a q × q orthogonal matrix, X is an n × n non-singular matrix, and Λ is a q × n diagonal matrix with diagonal elements The only elements of the m × n matrix ∆ that are possibly non-zero are One convenient property of the GSVD is that ∆ T ∆ + Λ T Λ = I n [3, p. 104]. The actual structure of ∆ depends upon the dimensions of A. If m ≥ n, then ∆ is a diagonal matrix. Let a k-diagonal matrix be defined as an m × n matrix with nonzero entries given by eq. (13) when m < n (which implies k = 0). In other words, a k-diagonal matrix is a rectangular matrix where the only entries that are possibly nonzero are located on the kth upper diagonal. With this definition, ∆ is a k-diagonal matrix if m < n.
An application of the decompositions in eq. (12) is that the normal equations corresponding to eq. (10) can be simplified as Hence, due to the invertibility of X.
The assumption that the block system matrix in eq. (10) has full column rank means that ∆ T ∆ + α 2 Λ T Λ is non-singular for α > 0. Then, the solution of eq. (14) is given by where Y is the inverse of X T . A representation of x(α), similar to eq. (5), is given by with δ = diag(∆ T ∆) and λ = diag(Λ T Λ) (the square roots being applied element-wise).

Discrete Cosine Transform
While the GSVD is useful for analyzing problems with a general matrix A, for practical image deblurring problems it is more efficient to use the 2D discrete cosine transform (DCT).
Assuming that reflective boundary conditions are applied, then both A and L have the same block structure and the DCT can be used to simultaneously diagonalize A and L. This special block structure is called BTTB + BTHB + BHTB + BHHB in [24], with the "T" and "H" standing for Toeplitz and Hankel, respectively. Theorem 2.1 describes how a simultaneous diagonalization of these matrices is related to their GSVD.
Theorem 2.1. [8] Given the simultaneous diagonalization of the n × n symmetric matrices A = C T ∆C and L = C T ΛC where C is orthogonal and null(A) ∩ null(L) = {0}, then the simultaneous decomposition is related to a GSVD which is given by A = U ∆X T and and P is the permutation matrix which reorders the entries in ∆S −1 in increasing order.
Proof. We begin by relating the decomposition of each matrix with a standard GSVD decomposition, eq. (12). First for A, we set C T ∆C = U ∆ X T , where U is orthogonal, X is invertible, and the entries in ∆ are in increasing order 0 ≤ ∆ 1,1 ≤ . . . ≤ ∆ n,n ≤ 1. Rearranging terms gives ∆ = U T C T ∆C X −T . Doing the same for L = V Λ X T , where the entries of Λ are in decreasing order, 1 ≥ Λ 1,1 ≥ . . . ≥ Λ n,n ≥ 0, we have Λ = V T C T ΛC X −T . Moreover, using the GSVD A T A + L T L = X( ∆ T ∆ + Λ T Λ) X T = X X T , while using the simultaneous decomposition we have A T A + L T L = C T (∆ T ∆ + Λ T Λ)C. Hence using the invertibility of C we have where S = C X. But now from the left hand side of the equation, S is symmetric and we have S 2 = ∆ T ∆ + Λ T Λ and since null(A) ∩ null(L) = {0}, the matrix on the right is diagonal with positive entries. Thus we can define the positive square root S = √ ∆ T ∆ + Λ T Λ, which is also diagonal with positive entries. Moreover we can associate the individual entries of the matrix ∆S −1 with cos(θ i ) for some i and the entries of ΛS −1 with sin(θ i ) for some i. Hence if we apply a permutation matrix P to each side of ∆S −1 so as to maintain diagonal structure, but to also reorder the entries from small to large, then the same similarity transformation applied to ΛS −1 will sort the entries in the opposite order, and we have a matrix of cosine entries and a matrix of sine entries, which we define to be ∆ = P T ∆S −1 P and Λ = P T ΛS −1 P , respectively. Clearly we have the standard normalization, ∆ T ∆ + Λ T Λ = I n . Moreover, we can take and we have the given GSVD for the matrix pair.
The condition null(A) ∩ null(L) = {0} in theorem 2.1 is necessary for the property ∆ T ∆ + Λ T Λ = I n . If null(A) ∩ null(L) = {0}, then the matrix on the left in ?? is singular and this contradicts the invertibility of C and Z on the right side. As an alternative, we could require that ∆ T ∆+Λ T Λ = I n where I n is a modified identity matrix that has some zero diagonal elements. This generalization, as well as a conversion involving Kronecker products, is described in Chapter 2 of [8]. Furthermore, equipped with theorem 2.1, it is clear that the terms in eq. (17) can be rewritten in terms of the DCT simultaneous decomposition.

Multi-parameter Tikhonov regularization
A more general approach to regularization replaces the single regularization parameter by a MP vector α = [α (1) , . . . , α (p) ] T . Here, we follow the approach in [10] by defining P vectors w (p) ∈ R n that contain non-negative weights which satisfy Defining windows W (p) = diag w (p) for p = 1, . . . , P , we have P p=1 W (p) = I n . Considering first the Tikhonov regularization with L = I n , a regularization parameter α (p) can be selected for each window W (p) so that a windowed regularized solution can be constructed as which corresponds to eq. (9) when P = 1. In this framework, the windows must be selected before choosing corresponding regularization parameters.
The condition given by eq. (21) means that for each j = 1, . . . , n, there is exactly one p ∈ {1, . . . , P } such that w (p) j = 1. When working with non-overlapping windows, the pigeonhole principle [15] can be used to show that there will exist indices p such that w (p) = 0 if P > n. Perhaps the simplest way of choosing the components of w (p) is to first choose P +1 partition values ω (0) ≥ . . . ≥ ω (P ) such that ω (0) ≥ s 1 and s n > ω (P ) , then for p = 1, . . . , P , we set There are some advantages of using eq. (22). One advantage is that singular values of similar magnitude are grouped together. Another advantage is that the MP estimation functions to be discussed in section 3 decouple into linear combinations of single parameter functions.
Choosing ω (1) , . . . , ω (P −1) to be the P − 1 linearly-spaced or logarithmically-spaced points between s 1 and s n and then setting ω (0) = s 1 and ω (P ) < s n is an example of how to use eq. (22). Partition values ω (0) ≥ . . . ≥ ω (P ) can also be used to generate overlapping windows. For example, cosine windows are defined in [10] by using midpoints The first and P th weight vectors can be defined by respectively. See [8] for a proof that such cosine windows satisfy eq. (19). For generalized Tikhonov regularization (L = I n ), a windowed regularized solution similar to eq. (20) can be obtained using the GSVD. Given the decompositions A = U ∆X T and L = V ΛX T described in section 2.2 and windows W (p) for p = 1, . . . , P , a corresponding windowed regularized solution is with Y again being the inverse of X T . In terms of the filter functions eq. (18), the windowed solution eq. (23) can be written as Care must be taken when selecting the weight vectors w p in light of the fact that the generalized singular values are arranged in ascending order. If we extend the notation Φ(α) and Ψ(α) introduced in section 2.2 to define the n × n diagonal matrices then eq. (23) can also be written as

Multiple data sets
We now consider the situation where we have a collection of data sets The vectors d (r) , b (r) , η (r) , and 0 (r) have length m r , while the vector x (r) has length n r . The system matrices A (r) and covariance matrices Σ (r) are thus m r ×n r and m r ×m r , respectively. We also assume that the random vectors {η (r) } R r=1 are mutually independent. For given regularization parameters α (r) and penalty matrices L (r) of dimension q r × n r , Tikhonov regularization can be performed to produce regularized solutions x(α (r) ) that minimize the functionals Typically, a parameter selection method is utilized to select the regularization parameter for each system. Instead, let d be the vector formed by vertically concatenating the data sets {d (r) } R r=1 and define the functional where α is a single regularization parameter. This construction implies that x, b, η are vertical concatenations of vectors in {x (r) } R r=1 , {b (r) } R r=1 and {η (r) } R r=1 , respectively. The matrices A and L are block diagonal matrices with diagonal blocks {A (r) } R r=1 and {L (r) } R r=1 , respectively. By the definition of the 2-norm and the construction of eq. (25), we can also write The advantage of regularizing via eq. (25) is that we only have to select one parameter instead of R parameters (one for each data set). Assumption 1 summarizes the set-up established in eq. (24).
) with the η (r) being mutually independent. The vectors b (r) , d (r) , and η (r) are of length m r and x (r) is of length n r .
To obtain an explicit representation of x( α), the GSVD can be utilized. A GSVD of each pair (A (r) , L (r) ) can obtained for r = 1, . . . , R so that Using the tilde notation introduced in section 1.1, we express A and L as Analogous to representation eq. (15) for a single system, the solution x( α) can be written as where Y is the inverse of X and d = U T d. The filter function notation from section 2.2 can be used to define , otherwise for r = 1, . . . , R, as well as the n r × n r diagonal matrices Concatenating Φ (r) (α) and Ψ (r) (α) via the block diagonal representation, eq. (26) becomes Some of the parameter estimation methods considered in this paper rely upon the statistical properties of the noise in the data, so the statistics of η will be addressed. Though eq. (24) indicates that the random vectors {η (r) } R r=1 are assumed to have zero mean, we can relax this assumption so that η (r) ∼ N (µ (r) , Σ (r) ) for all r = 1, . . . , R. The distribution of η is then given by lemma 2.1, which follows from the properties of the multivariate normal distribution.
We conclude section 2.5 by enumerating and discussing the additional underlying assumptions that will be utilized in section 3.
Assumption 3. For all r = 1, . . . , R, we assume that m r = m. In other words, we assume that the size of each data vector d (r) is the same.
Assumption 4. We assume that there exist matrices ∆ ∈ R m×n and Λ ∈ R q×n such that It should be noted that assumption 2 could be relaxed so that Σ (r) is any diagonal matrix D (r) , in which case a whitening transformation C (r) could applied so that ξ (r) = C (r) η (r) ∼ N (0 m×1 , I m ). For example, one could use the zero-phase component analysis (ZCA) transformation C (r) = Σ (r) −1/2 [5].
The strength of the assumptions presented increase in accordance with their numbering. Assumption 1 describes the most general situation involving MD considered in this paper, where the size of each problem and the covariance of each zero-centered multivariate Gaussian noise vector are potentially distinct. Assumption 2 then requires that the covariance matrices be constant diagonal matrices, though the sizes of the problems remain unrestricted. In contrast, assumption 3 specifies that each data set much be the same size m; however, the size of the solutions may still be district (n r can differ for different values of r). A consequence of assumption 4 is that the pairs (A (r) , L (r) ) all have the same generalized singular values. Though assumption 4 is strong, it is implied by assumption 5. Assumption 5 is the strongest, in that it not only implies assumptions 3 and 4 but demands the solution size be the same and the system and penalty matrices be the same for all problems. A useful consequence of assumption 5 is that Φ (r) (α) = Φ(α) and Ψ (r) (α) = Ψ(α) for all r = 1, . . . , R and α > 0.

Multi-parameter regularization for multiple data sets
We now consider the case where MP regularization is applied to multiple data sets (MDMP). Letting α (r) = [α (r,1) , α (r,2) , . . . , α (r,Pr) ] be the P r regularization parameters used for MP regularization applied to the rth system described by eq. (24), we can independently construct regularized solutions where W (r,p) is the pth window for the rth system. Each system can have its own set of windows, meaning that there are a total of R r=1 P r regularization parameters. The primary assumption we make moving forward, however, is that the number of windows W (r,p) for each system is the same, i.e P r = P for all r = 1, . . . , R. This implies that α (r) are all vectors of length P and there are a total of RP parameters. A stronger assumption would require windows are all the same, which is described by assumption 6 and used in section 3.
Analogous to the single parameter case for MD, we define x win ( α) as the vertical concatenation of the x (r) win ( α) for r = 1, . . . , R, where α ∈ R P is a single vector of parameters. We then have where W (p) = diag(W (1,p) , . . . , W (R,p) ) has R diagonal blocks. As the final extension of notation from section 2.5, we define for r = 1, . . . , R and p = 1, . . . , P . With the diagonal matrices the appropriate block matrices can be formed: Using Φ α (p) for the pth window, representation eq. (27) can be written as

Parameter selection methods
As stated in section 1, there are a variety of methods for selecting regularization parameters and we illustrate some representative examples. The UPRE and MDP methods are popular choices for situations in which the variance of the noise in the data is known. In contrast, the GCV and L-curve methods are common choices when the variance is unknown. However, all four methods involve the norm of the regularized residual r(α) [50, p. 98]. For the standard single parameter, single data case, the regularized residual is defined as The UPRE and GCV methods require the trace of the influence matrix A(α); the single parameter, single data set version is given by If the GSVD is used with the filter factor matrix Φ(α), we can express A(α) as where the last equality follows from the identity ∆ T = ∆ T ∆∆ † [8].
To obtain the functions associated with the UPRE, MDP and GCV methods, we first discuss how to obtain expressions for the regularized residual and trace of the influence matrix, under the MP and MDMP situations in section 3.1. We then give in sections 3.2 to 3.4 the explicit functions that are obtained in the context of the UPRE, MDP and GCV methods.

Regularized Residual and Traces for MDMP Solutions
To use the regularization methods we need to obtain expressions for the regularized residual and influence matrix for the MDMP case. We extend the notation introduced in section 2.4 for MP regularization involving a single data point to define a windowed regularized residual r win (α) = Ax win (α) − d where the windowed regularized x win (α) is given by eq. (23) assuming the use of the GSVD. Similarly, we generalize the windowed regularization matrix defined in [10] so that for windowed spectral regularization, the influence matrix A win (α) is If instead we have MD but each is used with a scalar α, we can extend the notation from section 2.5 so as to define a regularization solution and influence matrix A GSVD representation of x( α) is provided by eq. (26) and the diagonal blocks of A( α) are For the most general MDMP case, we combine the notations of sections 2.4 to 2.6 to construct a regularized residual x win ( α) can be represented by eq. (27), while the diagonal blocks of A win ( α) are Both W (r,p) and Φ (r) ( α (p) ) are n r × n r matrices.
We will now present representations of the necessary norm and trace terms that are used in discussions of the UPRE, MDP, and GCV methods adapted for this case. These representations are presented for MP regularization of MD since this is the most general situation. The first representation is given in theorem 3.1, which concerns the norm of the regularized residual. with k r = n r − m r .
Proof. We first let M = R r=1 m r , which is the length of r win ( α). Substituting eq. (29) into the definition of r win ( α), we have Using the 2-norm and the block structure of the matrices, we can then write We must now consider the two cases for each r = 1, . . . , R. If m r ≤ n r , then let k r = n r − m r and so that If m r > n r instead, then let k r = m r − n r . The matrix within the norm has the block form Thus, the norm becomes As a function of α ∈ R P + , the limiting behavior of r win (α) 2 2 provides insight into how best to deal with minimizing functions that involve r win (α) 2 2 . Lemma 3.1 shows that r win (α) 2 2 can be bounded above by the (scaled) norm of the data. Lemma 3.1. Given P ∈ {1, . . . , N } windows, Proof. From theorem 3.1 with R = 1 and A ∈ R m×n , we have that , m > n.
with k = n − m. Without loss of generality we assume that m ≤ n, so From eq. (28), we know that 0 ≤ ψ j+k (α (p) ) ≤ 1 for all α (p) ∈ R + . Additionally, 0 ≤ w (p) j+k ≤ 1 for all p = 1, . . . , P . Therefore we have with the last equality following from d = U T d with orthogonal U .
In practice, lemma 3.1 does not provide a tight bound since many of the terms in the sum over p are less than one. If non-overlapping windows are used, then a tighter bound could be obtained by noting that many of the summands are identically zero. However, lemma 3.1 does provide some insight into the numerical behavior of the functions using r win (α) 2 2 , which will be discussed in section 4.
We can also develop an analog of the windowed regularized residual as applied to the average d = 1 R R r=1 d (r) under assumptions (5) to (6). Defining the windowed regularized residual applied to d as where the windowed regularized solution x win ( α) is with d = U T d, the following corollary then applies.
Corollary 3.1.1. Under assumptions (1) to (6), for all α ∈ R P + we have that Proof. Under assumption 5, m r = m, n r = n, and (A (r) , L (r) ) = (A, L) for all r = 1, . . . , R. Without loss of generality, suppose that m ≤ n and let k = n − m. Then using eq. (33) and eq. (34), we have We also have Thus, through a change of summation where the last equality follows from theorem 3.1 in conjunction with assumption 6.
Proposition 3.1 describes how r win ( α) 2 2 can be decomposed when working with nonoverlapping windows.  where r (r,p) with k r = n r − m r .
Proof. Given r ∈ {1, . . . , R}, we assume without loss of generality that m r ≤ n r . Using theorem 3.1, we can write r with k r = n r − m r . Since the weight vectors {w (r,p) } P p=1 satisfy eq. (21), for each index j there exists exactly one index p such that w (r,p) j = 0. Thus, the sum over p has only one nonzero summand for each index j, meaning we can write Therefore r Proposition 3.1 means that if regularization is being performed with non-overlapping windows, then the norm of the regularized residual can be written as a sum of norms of residuals specific to each window. Now, theorem 3.2 provides a representation of the trace term that is needed in section 3.
Proof. The diagonal block structures of both Σ = diag Σ (1) , . . . , Σ (R) and A win ( α) allow the trace to be written as a sum of traces: From assumption 2, we have that Σ (r) = σ 2 r I mr for all r = 1, . . . , R. Furthermore, representation eq. (32) of the diagonal blocks of A win ( α) and the similarity invariance of the trace operation allow us to write Using k r = min{0, n r − m r }, we have Therefore, With the inclusion of assumption 5 and assumption 6, we can make a statement regarding the traces of the covariance and influence matrices similar to that of corollary 3.1.1. We first define Σ as the covariance matrix of the averaged noise η = 1 R R r=1 η (r) for the averaged data d. Since the random vectors {η (r) } R r=1 are mutually independent and have mean zero, E(η) = 0 m×1 and the evaluation of Σ is reduced to  (1) to (6), we have that Proof. We prove item (i) only, since item (ii) follows directly from eq. (35) and the block structure of Σ. Using eq. (35), for all α ∈ R P + we have trace Σ (r) A win (α) .

Assumption 2 then gives
Shifting to 1 R 2 trace Σ A win (α) , assumptions 5 and 6 imply A In contrast to proposition 3.1, it is not necessary to have non-overlapping windows in order to decompose trace terms into traces for each window. The window decomposition of the trace is described in proposition 3.2.
Proof. From theorem 3.2, we can write trace Σ (r) A (r) with k r = min{0, n r − m r } for each r = 1, . . . , R. Changing the order of summation then gives where A (r,p)

Unbiased Predictive Risk Estimator
The UPRE method was developed in 1973 by Mallows and considers the statistical relationship between the regularized residual eq. (30) and the predictive error p(α) = A(x(α) − x).
We first present the UPRE function for Tikhonov regularization under the more general condition that η ∼ N (0, Σ). Analogous to the standard UPRE function, we can redefine F UPRE (α) as See [8] for a statistical derivation of eq. (37). The standard UPRE function eq. (36) is recovered from eq. (37) if Σ = σ 2 I m (which is assumption 2 with R = 1).
where F (r) UPRE win Proof. We first let M = R r=1 m r , which is equal to mR by assumption 3. Lemma 2.1 implies η ∼ N ( 0, Σ) with 0 ∈ R M , and so the more general UPRE function eq. (37) can be directly applied to produce As used in the proofs of theorems 3.1 to 3.2, the structure of r win ( α), Σ, and A win ( α) allow for eq. (40) to be written in terms of sums over r: .
Then for all α ∈ R P + , The proof of proposition 3.4 proceeds by using M = mR with corollary 3.

Morozov's discrepancy principle
Similar to the UPRE method, the MDP method relies on knowledge of the variance of η. If η ∼ N (0, σ 2 I m ), then the MDP function is The function eq. (42) stems from the observation that if The MDP method selects α MDP as the zero of F MDP (α). Before extending the MDP method to account for MD, some comments on the behavior of the MDP function eq. (42) are warranted. The term 1 m r(α) 2 2 is a monotone increasing function of α for α > 0 [50]. The monotonicity of 1 m r(α) 2 2 does not guarantee, however, the existence of a zero of F MDP (α). If the selected value of σ 2 is too large, then it is possible that F MDP (α) < 0 for all α > 0 and a root will not exist. This can be attributed to the limiting behavior of 1 m r(α) 2 2 , which is described in [8]. If the selected value of σ 2 is larger than 1 m d 2 2 , F MDP (α) will not have a root for α > 0 and the MDP method fails to select a regularization parameter. Sometimes a safety parameter > 0 is introduced to modify the MDP function to to account for root-finding difficulties [3,18], though selecting an appropriate value of is an ad hoc process and depends on the confidence of σ 2 as the true noise. The original MDP function is recovered from eq. (43) when = 1. The MDP function for the more general assumption that η ∼ N (0, Σ) will now be presented. Since E( 1 m r(α) 2 2 ) = E( 1 m η 2 2 ) and E( η 2 2 ) = trace(Σ), the MDP function is redefined to be Equation (44) can be applied to data {d (r) } R r=1 .
Proof. Using lemma 2.1 for η with M = R r=1 m r = mR, the MDP function eq. (44) can be applied to define By defining the individual MDP functions by eq. (45), the MDP function eq. (46) for the large system can be written as an average by exploiting the block structure of A(α) and Σ: The MDMP MDP method then defines α MDP as the zero of F MDP win ( α). Analogous to the UPRE method, eq. (45) is equivalent to eq. (42) under assumption 2. A safety parameter can also be included in the trace terms of eq. (46) and eq. (45) for more control over selected parameters.
As with the MDMP UPRE method, the MDMP MDP method is distinct from the MDP method as applied to the averaged data d. Proposition 3.6. Under assumptions (1) to (6), let Then for all α ∈ R P + , The proof of proposition 3.6 is nearly identical to that of proposition 3.4.

Generalized cross validation
Unlike the UPRE and MDP methods, the GCV method does not depend upon knowledge of the noise variance. For a single data set, the GCV function is defined as The GCV method determines α such that α GCV = arg min α>0 F GCV (α). In contrast to the UPRE and MDP methods, the GCV function for windowed regularization of MD is not equal to the mean of individual GCV functions when only assumptions (1) where F (r) GCV win ( α) is defined by eq. (48).
Proof. Again we let M = R r=1 . Assumption 5 implies that A (r) win ( α) = A win ( α) for each r = 1, . . . , R, as well as M = mR. Having the same influence matrix for each r is crucial because of the trace term contained in the denominator of the MDMP GCV function From theorems 3.1 to 3.2, we can then write The MDMP GCV method defines α GCV = arg min α>0 F GCV win ( α). Note that eq. (48) is equivalent to eq. (47) without assumption 2, which is to be expected since the GCV method does not rely on knowledge of the noise variance.
A result analogous to proposition 3.4 and proposition 3.6 can be obtained for the MDMP GCV function using corollary 3.1.1 and corollary 3.2.1. Then for all α ∈ R P + , As a closing remark regarding the MDMP GCV method, it is important to note that assumption 5 was needed for result proposition 3.7. In contrast, no such assumption is necessary for the corresponding MDMP UPRE and MDMP MDP results, which respectively are proposition 3.3 and proposition 3.5.

Validation
To evaluate the effectiveness of both the MD and MDMP parameter selection methods described in section 3, two test problems were considered. For both problems, data vectors were considered as a set of training data and resulting parameters were then applied to validation data set. A 1D problem serves as a proof of concept for the MD methods and uses MRI data that is built into MATLAB ® , which is shown in fig. 1. The second problem  Figure 1: MRI data formed by reformatting the MATLAB ® built-in MRI data. Each vector was blurred using a Gaussian kernel (centered at zero with variance 1 (each pixel representing a unit square in the continuous problem)). Normal noise was added to produce a specific SNR (eq. (50)); the resulting vectors are shown at the bottom. The dimension of all three images is 256 × 536. The SNR of each column was selected from a uniform distribution between 6 and 7.
is two-dimensional (2D) and utilizes the images in fig. 2 of the planet Mercury obtained by the MESSENGER space probe 1 . For both problems, the signal-to-noise ratio (SNR) is used as a measurement for noise content in the images and is given by SNR = 10 log 10 P signal P noise .
In the discrete setting, the average power P of a vector x of length n is defined as 1 m x 2 2 . Using this definition for vectors b and η, P signal = 1 m b 2 2 and P noise = 1 m η 2 2 and so the quotient in the logarithm is b 2 2 / η 2 2 . If b is a matrix representing an image, in which case η is a realization of a random matrix, the 2−norms can be replaced by the Frobenius norm.
For a basis of comparison, parameters were also selected as minimizers of the learning function where F (r) Best win Note that this definition requires that the true solutions, {x (r) } R r=1 , are known for generating the trained data. Regularization parameters chosen as minimizers of eq. (51) are "best" in the sense of minimizing the mean squared errors of the regularized solutions x (r) win (α); the use of eq. (51) was considered in [11]. One could also find minimizers α (r) of eq. (52) for each r = 1, . . . , R, which would produce parameters that are "best" for their own data set. In the results we use Best to indicate results that are found using the learning function eq. (52). Moreover, in all figures and tables we use the shorter versions UPRE, MDP and GCV,for the methods that are implemented for MD or MDMP, with an explanation in the captions, rather than presenting cluttered legends.

One-dimensional problem
The 1D test problem uses the MRI data built into MATLAB ® , which can be accessed by the command load mri.mat. Five horizontal slices were selected and reformatted as a single image ( fig. 1) of dimensions 256 × 536. The default dimension of each horizontal MRI slice accessible using load mri.mat is 128 × 128. Linear interpolation was used to double the number of rows; the number of columns of the concatenated MRI slices was trimmed to to eliminate leading and trailing zero columns. The columns of the image were then multiplied by a symmetric Toeplitz matrix, approximating a Fredholm integral equation of the first kind with a zero centered Gaussian kernel of variance 1 (each pixel representing a unit square in the continuous problem); explicitly, the continuous kernel is k(y) = exp (−y 2 /2). The resulting MRI image is vertically blurred; see fig. 1. SNR values were selected randomly between 6 and 7 for each blurred column, and realizations of normal noise with the corresponding variances were added to blurred columns. The first 536/2 = 268 columns of fig. 1 served as the training set, while the remaining 268 columns serve as a validation set. For the penalty matrix L, the standard one dimensional first order derivative matrix of size 255 × 256 was used. It is everywhere 0 except for a 1 on each super-diagonal entry and −1 for each diagonal entry. The resulting system matrix in eq. (10) has full column rank, and so applying the normal equations in terms of the GSVD results in unique solutions. Figure 3 illustrates a comparison of the methods in terms of the regularization parameters α selected as the number of training vectors increases. When the number of training vectors is small (e.g. between 1 and 10 from fig. 3), the parameters determined by all four methods change significantly. This can be attributed to the fact that all of the methods find a parameter that is either a root or a minimum of an average of functions. For a small number of training vectors, each additional vector has more influence over the behavior of the MD function. The parameters stabilize as the number of training vectors reaches a certain point; in the problem being considered here, the parameters stabilize by about 40 training vectors.
Before looking at the relative errors of the regularized solutions, another observation regarding fig. 3 can be made. While the parameters determined using the learning function eq. (51) and the MD UPRE and MD GCV methods are close (the maximum difference in the parameters does not exceed 0.5), the MD MDP method with = 1 consistently selects significantly larger parameters. A direct consequence is that the regularized solutions have higher relative errors. Figure 4 shows the relative errors of the regularized solutions corresponding to the parameters from each method. Instead of using a relatively small number of training vectors, such as 40 out of 268 in fig. 3, the full training set was used to determine a parameter by each method. These four parameters were then used to generate a regularized solution for each vector in the validation set, and the relative errors were computed against the true solutions. First, the relative errors obtained by the learning and MD UPRE and MD GCV methods are quite similar; the means of the relative errors are close, and there is a collection of upper outliers. In contrast, the MD MDP method has a slightly larger mean relative error. The most striking visual characteristic from fig. 4, however, is the number of outlier relative errors generated by the MD MDP. The largest relative error obtained using the MD MDP method is about twice the largest relative error of the other three methods. The 1D test problem demonstrates that the MD methods have potential for selecting a viable regularization parameters that can be applied for multiple sets of data. The MD UPRE and MD GCV methods performed competitively against the learning method, which relies on knowledge of the true training solutions. As stated previously, the success of the MD methods is predicated on the similarity of the data sets being considered. This is not an unreasonable condition because data from a given experiment would hopefully vary little without significant change in the experimental set-up.

Two-dimensional problems
The data sets of the 2D test problem consist of images of size 256 × 256. A total of 16 images were used and split into training and validation sets containing 8 images each. To obtain the 16 images from the 512 × 512 Mercury images in fig. 2, 8 images were randomly chosen and two 256 × 256 subimages of each image were randomly selected. The 16 images were then split into training and validation sets. Another validation set of images, shown in fig. 5, was used that consisted of built-in MATLAB ® images.
A 256 × 256 point spread function (PSF) was formed using a discretization of the zero centered, circularly symmetric Gaussian kernel The parameter ν controls the width of the Gaussian kernel. Choosing k(x, y) to be circularly symmetric is for convenience; a Gaussian kernel with different width parameters for the x and y directions can still be used to construct a PSF that is doubly symmetric for diagonalization via the DCT [24]. In regards to the value of ν, values ν = 4, 16, and 36 correspond to blurring that is referred to as "mild," "medium," and "severe" in [18], respectively. The corresponding PSFs were discretely convolved with each image as a means of blurring. SNR values of 10 and 25 were used to construct realizations of normal noise that were added to the blurred images to create the data. For one choice of the penalty matrix L, we used the appropriately structured version of the discrete negative Laplacian operator, which is an approximation of the continuous Laplacian operator [14,31]. For the second penalty matrix we used L = I. The structure of A and L allows for simultaneous diagonalization using the DCT for numerical efficiency (see section 2.3).
In regards to the spectral windowing, at most two windows were used (corresponding to the use of just two parameters in the MP estimators). Two non-overlapping logarithmic windows and two overlapping logarithmic cosine windows were considered; the decision to use logarithmic spacing is supported by how the ordered spectral components decay, which is shown in fig. 6. The use of three or more windows was not considered for the sake of computational savings and because the use of more than two windows was not advocated in [10]. The use of three windows did not produce noticeable benefits in our test problems as well. Lastly, results for the MDMP MDP method as applied to the 2D problem have been omitted since results depend upon the tuning of the safety parameter .
For all of the numerical tests involving two parameters, the relationship between the parameters mostly depends upon the choice of penalty matrix. For L = I, the first parameter α (1) dominates the resulting regularized solutions even though the magnitude of α (1) remains small (less than 0.2 in fig. 7). The second parameter α (2) tends to be a few orders of magnitude larger than α (1) but is localized; fig. 7 shows α (2) bounded between 9.9 and 10.4 for all three methods. Figure 8 shows the contrast in parameter behavior when the L is chosen to be the Laplacian. The magnitude of α (1) increases while the magnitude of α (2) decreases. For all three methods, the range of both parameters also increases. The relative errors, however, of the associated regularized solutions are consistent across methods, for each choice of L.
The use of overlapping or non-overlapping windows influences the degree of interdependence between the two parameters. Figures 8 to 9 present the results of using overlapping and non-overlapping logarithmic windows. When using overlapping windows ( fig. 8), the ranges of both parameters are smaller than those for non-overlapping windows ( fig. 9). The exception is the behavior of α (2) exhibited by the use of F Best win (α) in fig. 9, which shows all of the second parameters grouping near 10. This grouping of the second parameter was   Figure 6: Decay of the ordered operator spectra in relation to width of the Gaussian kernel and choice of penalty function. The vertical lines indicate the partition points when using two (linear or logarithmic) windows. Based on the decay of the spectra, logarithmic rather than linear windows were selected for the numerical tests. also exhibited in other numerical tests involving non-overlapping windows. The grouping behavior can be explained by the choice of an upper bound during the minimization process; in the case of fig. 9, the upper bound was chosen near 10. The calculated gradient of the F Best win (α) is too small to resolve a minimum in the direction of α (2) and thus the minimization process determines the minimizers near the specified boundary.

Decay of spectral components for choices of blur amount and penalty matrix
In regards to the MDMP methods, which select P parameters using R data sets, the trend of the parameters as R increases is similar to that displayed by the 1D case ( fig. 3) for the first parameter only. Figure 11 shows that α (1) stabilizes using four data sets or more for all three methods. While the MDMP UPRE and MDMP GCV methods produce stable values of α (1) near the same value of R, the learning method produces values that oscillate between the specified upper bound (20 in this case) and values similar to the other two methods. This failure to stabilize can most likely be attributed to the difficulties described earlier regarding minimization in the direction of α (1) . Since the value of the second parameter often has little effect on the evaluations of the functions, the oscillatory nature of α (1) does not negatively effect the relative errors of the regularized solutions. Table 1 details the mean percent relative errors of solutions obtain using parameters from each MDMP method. Even for the limit number of training sets (2 through 8), the errors decrease as R increases.
It is interesting to note that the average relative errors of solutions obtained for parameters applied to the second validation set ( fig. 5) were less than those of either the training or first validation set. The similarity of errors between the training and first validation set is to be expected since both sets consist of images from the same parent set, which are subimages  Figure 7: Boxplots of the parameters α (1) and α (2) , as well as the relative errors of the corresponding regularized solutions. The MDMP UPRE and MDMP GCV methods select similar values of α (1) ; all three methods select similar values for α (2) . The distributions of the relative errors are similar for all three methods, with MDMP UPRE and MDMP GCV methods producing results slightly greater than obtained using the learning function.   Figure 9: Boxplots of the parameters α (1) and α (2) , as well as the relative errors of the corresponding regularized solutions. By choosing non-overlapping logarithmic windows, the learning function selects α (2) at the upper limit set during the minimization process while the MDMP UPRE and MDMP GCV methods provide values having a wider distribution. selected from fig. 2. The superior (reduced) errors calculated for the second validation set is consisted throughout most numerical configurations, though the difference in relative errors between the use of the different validation sets never exceeded ten percent. Additionally, some solutions from the second validation set have relative errors greater than those in either the training or first validation set, but these are offset by smaller errors during averaging. Figure 10 presents two examples of images from the second validation set that have differing relative errors.  Figure 10: Two samples from the second validation set, with ν = 36, an SNR of 10, two log cosine windows and the Laplacian penalty matrix. From left to right for each sample are the true solution, the blurred image, the blurred image after noise was added, and the regularized solution. The two parameters were selected using the MDMP UPRE method with R = 8 (the entire training set).
Finally, taken together, figs. 3 to 11 suggest that the obtained parameters stabilize fairly fast in relation to the number of data sets used for parameter estimation, though there can be issues with minimization in the direction of additional parameters.

Conclusions and future work
We have shown that the UPRE, MDP, and GCV methods can be extended to accommodate regularization parameter estimation using multiple data sets and for both single and multiple parameters, for generalized Tikhonov regularization. The UPRE and MDP are representative estimators that assume the knowledge of the variance of mean zero Gaussian noise in the data,  Figure 11: MDMP parameters obtained from using two log cosine windows as the number of data sets increases from 1 to 16. For all three methods, α (1) stabilizes while α (2) stabilizes only for the MDMP UPRE and MDMP GCV methods; the parameters from the learning method oscillates between the specified upper bound and parameter similar to the other two methods.
while no additional assumptions are required for the GCV estimator. The most general forms of functions associate with these methods are eq. (40), eq. (46), and eq. (49), respectively. While the corresponding functions for MDMP UPRE and MDMP MDP methods can be written as an average of the individual functions associate with each data set, this is not possible for the MDMP GCV method without stronger assumptions (see section 3.4). None of the three MDMP methods require knowledge of true solutions unlike the learning approach defined by eq. (51). The presented numerical experiments for 1D and 2D signal restoration demonstrate that the MDMP methods can perform competitively with the learning approach that requires knowledge of true signals for training the data. Further, it is also demonstrated that the parameters obtained from a specific training set of validation images can also be used for a set of different testing images, provided that the general noise characteristics are the same. The approach extends immediately for any estimator which relies only on an approximation for the regularized residual and the trace of the influence matrix, and the general idea can be modified to address estimators requiring other terms, such as an augmented regularized residual use for the χ 2 estimator described in [34,35]. Although the derivations are presented for the case in which there is a known mutual decomposition of the model and penalty matrices (A and L), because it is only needed to obtain estimates of the required terms, the approach can be extended for any iterative method which yields suitable estimates, e.g. [12,42].
Further directions for future investigation include (i) determining a priori the necessary number of data sets for a stabilization of selected parameters, and (ii) investigating the implementation of MDMP methods in the framework of more general regularizers, such as total variation regularization which is imposed to preserve sharp edges in images. Though the choice of the identity and the Laplacian as penalty matrices are not conducive to sharp regularized images, their use here lays the foundation for more intricate regularization approaches, as applied for TV regularization in the context of iteratively reweighted and augmented Lagrangian algorithms.