# A Regularized Alternating Least-Squares Method for Minimizing a Sum of Squared Euclidean Norms with Rank Constraint

Authors: Pablo Soto-Quiros
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

Minimizing a sum of Euclidean norms (MSEN) is a classic minimization problem widely used in several applications, including the determination of single and multifacility locations. The objective of the MSEN problem is to find a vectorx such that it minimizes a sum of Euclidean norms of systems of equations. In this paper, we propose a modification of the MSEN problem, which we call the problem of minimizing a sum of squared Euclidean norms with rank constraint, or simply the MSSEN-RC problem. The objective of the MSSEN-RC problem is to obtain a vector x and rank-constrained matrices A1,⋯,Ap such that they minimize a sum of squared Euclidean norms of systems of equations. Additionally, we present an algorithm based on the regularized alternating least-squares (RALS) method for solving the MSSEN-RC problem. We show that given the existence of critical points of the alternating least-squares method, the limit points of the converging sequences of the RALS are the critical points of the objective function. Finally, we show numerical experiments that demonstrate the efficiency of the RALS method.


---

## Body

## 1. Introduction

The problem of minimizing a sum of Euclidean norms (MSEN) is an optimization problem whose goal is to find a vectorx∈ℝn such that it minimizes a sum of Euclidean norms of systems of equations Ajx=bj, where Aj∈ℝm×n and bj∈ℝn, for j=1,⋯,p. Vectors bj represent the coordinates of existing facilities in the plane, and x represents the coordinates of a new facility in the position, which minimizes the weighted sum of distances from the new facility to each existing facility facilities ones. The multifacility location problem arises when the number of new facilities is more than one. The weighted sum of norms to be minimized then involves the distances between each pair of new facilities and the distances between the new and existing facilities.The MSEN problem was studied by Fermat in the seventeenth century, when eachAj is the identity matrix, m=n=2 and j=1,2,3. Here, the goal of this problem is to find the point in ℝ2 such that it minimizes the sum of distances from it to three given points. Moreover, various names are associated with the MSEN problem, for example, the general Fermat problem, the Weber problem, and the weighted single-facility location problem.The MSEN problem arises in many applications, such as the superresolution mapping based on spatial-spectral correlation to overcome the influence of linear and nonlinear imaging conditions [1]. Specifically, the spatial correlation is obtained using the mixed spatial attraction model based on the linear Euclidean distance. Other applications of the MSEN problem are the Euclidean facilities location problem applied in transportation and logistics and the Steiner minimal tree problem under a given topology [2–4].Several techniques solve the MSEN problem. For example, a predictor-corrector algorithm developed by Andersen et al. in [2] is derived from a primal-dual interior-point algorithm by applying Newton’s method directly to a system of nonlinear equations representing primal and dual feasibility and a perturbed complementary condition. Furthermore, the primal-dual algorithm developed by Qi et al. in [5] presents an augmented smoothing algorithm for solving a nonsmooth equation. By transforming the MSEN problem into a classic convex programming problem in conic form, Xue and Ye show in [6] that an ε-optimal solution of the MSEN problem can be estimated efficiently using interior-point algorithms.The MSEN problem is convex but not differentiable [3]. Moreover, it may happen that the accuracy of this problem is still not satisfactory because vector solution x is not a solution of all linear systems Ajx=bj, producing a higher value in the objective function of the MSEN problem.In this paper, we wish to develop and justify a new optimization problem that provides better-associated accuracy than the MSEN problem. We propose a new optimization problem based on the classical MSEN problem. This optimization problem is called the problem of minimizing a sum of squared Euclidean norms with rank constraint, or simply the MSSEN-RC problem. The goal of the MSSEN-RC problem is to obtain a vectorx and rank-constrained matrices A1,⋯,Ap such that they minimize a sum of squared Euclidean norms of p systems of equations Ajx=bj.The proposed MSSEN-RC problem is crucial because it extends the use of linear formsFjx=Ajx in the MSEN problem to bilinear forms Gjx,Aj=Ajx. The use of bilinear forms increases the accuracy of the objective function by increasing the number of parameters to be optimized. Augmenting the number of optimized parameters to improve accuracy is used in various minimization problems (see, e.g., [7–10]). Moreover, the objective function on the MSSEN-RC problem is differentiable, unlike the MSEN problem. This feature is essential because we can use a gradient method to estimate a solution and thus improve the prediction ability by updating the guess in the direction of your gradient, for example, the BFGS algorithm [11].There are two significant differences between the MSEN and the MSSEN-RC problems. First, the MSSEN-RC problem finds the optimal vectorx and rank-constrained matrices Aj, while the MSEN problem finds the optimal vector x. This modification is essential to obtain p equivalent linear systems Ajx=bj. Second, the MSSEN-RC problem considers a least-squares problem instead of a least-absolute problem as in the MSEN problem. This adaptation permits us to obtain a stable solution because the objective function in the MSSEN-RC problem is a differentiable function, unlike the objective function in the MSEN problem.Additionally, in this paper, we propose an algorithm to obtain a solution for the MSSEN-RC problem. This algorithm is based on the regularized alternating least-squares (RALS) method (In [12, 13], the RALS method is known as a proximal point modification of the Gauss-Seidel method) [13–15]. The RALS method is an iterative method that uses a regularization term in the classical alternating least-squares method [14]. This additional term does not address the degeneracy problem, however. Moreover, the limit points of the RALS method will be shown as the critical points of the original optimization problem and not of the regularized version (see Section 4.4). The RALS method is utilized in related fields, such as remote sensing and computer vision. For example, this iterative method is used to estimate the solution of a nonnegative coupled hyperspectral image superresolution model for image reconstruction [16].The RALS method is based on two known techniques. First, it uses an update to the classical Broyden Fletcher Goldfarb Shanno (BFGS) method as presented by Li and Fukushima in [17]. The BFGS method estimates a vector x that minimizes objective function f when matrices A1,⋯,Ap are known. The BFGS method is part of the Quasi-Newton family and is an iterative method for solving unconstrained nonlinear optimization problems (see Section 4.1). Second, the RALS method uses a solution to the generalized low-rank matrix approximation problem [18] to compute rank-constrained matrices A1,⋯,Ap that minimize f when x is known (see Section 4.2). We show that the limit points of the converging sequences of the RALS method are the critical points of f. Finally, in this paper, we present several numerical simulations that illustrate the advantages of the proposed RALS method to estimate a solution to the MSSEN-RC problem.The main motivation of this paper is to expand the theoretical framework developed around the classical MSEN problem (see, e.g., [2, 4]). Further, the main contributions of this paper are as follows:
(i)
We present a new optimization problem, so-called the MSSEN-RC problem. The MSSEN-RC problem generalizes the MSEN problem(ii)
A brief theoretical study on the MSSEN-RC problem is provided in this work(iii)
We develop an algorithm to estimate a solution for the MSSEN-RC problem. Additionally, we present its respective convergence analysis(iv)
Numerical experiments have been performed to demonstrate the effectiveness of the proposed algorithmThe paper is organized as follows. In Section2, we present the notation used in this paper. Additionally, we introduce the mathematical formulation of the MSEN problem. A theoretical study for the MSSEN-RC problem is developed in Section 3. Section 4 presents the RALS method to estimate a solution to the MSSEN-RC problem. Here, we explain in detail the two methods on which the RALS method is based: the BFGS method (in Section 4.1) and the low-rank matrix approximation technique (in Section 4.2). In Section 4.3, we analyze two approaches to compute regularization parameters. The convergence analysis of the proposed method is presented in Section 4.4. Section 5 presents numerical experiments based on random matrices. Finally, Section 6 contains some concluding remarks.


## 2. Preliminaries

### 2.1. Notations

Throughout this paper, we use the following notation. Letℝrm×n be the set of all real m×n matrices of rank at most r≤minm,n. Let ℝm×np be the p-ary Cartesian product over the set of all real m×n matrices, i.e.,
(1)ℝm×np≔X=X1,⋯,XpXj∈ℝm×n,∀j=1,2,⋯,p.We defineFr1,⋯,rpp⊆ℝm×np as the p-ary Cartesian product over the set of all real m×n matrices of rank at most rj≤minm,n, for j=1,2,⋯,p, i.e.,
(2)Fr1,⋯,rpp≔X=X1,⋯,XpXj∈ℝrjm×n,∀j=1,2,⋯,p.To simplify notation, we will omit the subscript inFr1,⋯,rpp, i.e., Fp=Fr1,⋯,rpp. We consider function f:ℝn×ℝm×np⟶0,∞ such that
(3)fx,A=∑j=1pfjx,Aj,where A=A1,⋯,Ap∈ℝm×np and fjx,Aj=Ajx−bj22. Based on equation (119) in [19], f is differentiable. The gradient vector ∇f:ℝn×ℝm×np⟶ℝn×ℝm×np is defined by
(4)∇fx,A=∂f∂xx,A,∂f∂A1x,A,⋯,∂f∂Apx,A,where
(5)∂f∂xx,A=2∑j=1pAjTAjx−bj,∂f∂Ajx,A=2Ajx−bjxT.Finally, it follows from Section 1.9 in [14] that a norm of ∇fx,A∈ℝn×ℝm×np is defined by
(6)∇fx,A=2∑j=1pAjTAjx−bj22+∑j=1p2Ajx−bjxTfr2,where ∥·∥fr is the Frobenius norm.


### 2.2. The MSEN Problem

Given matricesAj∈ℝm×n and vectors bj∈ℝm, for j=1,⋯,p, the MSEN problem finds a vector x^∈ℝn such that
(7)∑j=1pAjx^−bj2≔minx∑j=1pAjx−bj2.The MSEN problem emerges in many applications, such as the superresolution mapping, VLSL design, transportation and logistics, and the Steiner minimal tree problem. Furthermore, there are several algorithms to estimate optimalx^ in the MSEN problem, for example, the predictor-corrector method developed in [2], the primal-dual algorithm created in [5], and the interior-point technique presented in [6].Note that ifx^ is a solution of problem (7) such that
(8)∑j=1pAjx^−bj2=0,then x^ is a solution of p equivalent linear systems Ajx=bj. However, a solution to the problem (7) does not always satisfy the condition (8) because the accuracy of the MSEN problem is not yet satisfactory. To increase this accuracy, we also need to increase the number of parameters to optimize. Therefore, in the below section, we proposed a new minimization method called the problem of minimizing a sum of squared Euclidean norms with rank constrain, or simply the MSSEN-RC problem.


## 2.1. Notations

Throughout this paper, we use the following notation. Letℝrm×n be the set of all real m×n matrices of rank at most r≤minm,n. Let ℝm×np be the p-ary Cartesian product over the set of all real m×n matrices, i.e.,
(1)ℝm×np≔X=X1,⋯,XpXj∈ℝm×n,∀j=1,2,⋯,p.We defineFr1,⋯,rpp⊆ℝm×np as the p-ary Cartesian product over the set of all real m×n matrices of rank at most rj≤minm,n, for j=1,2,⋯,p, i.e.,
(2)Fr1,⋯,rpp≔X=X1,⋯,XpXj∈ℝrjm×n,∀j=1,2,⋯,p.To simplify notation, we will omit the subscript inFr1,⋯,rpp, i.e., Fp=Fr1,⋯,rpp. We consider function f:ℝn×ℝm×np⟶0,∞ such that
(3)fx,A=∑j=1pfjx,Aj,where A=A1,⋯,Ap∈ℝm×np and fjx,Aj=Ajx−bj22. Based on equation (119) in [19], f is differentiable. The gradient vector ∇f:ℝn×ℝm×np⟶ℝn×ℝm×np is defined by
(4)∇fx,A=∂f∂xx,A,∂f∂A1x,A,⋯,∂f∂Apx,A,where
(5)∂f∂xx,A=2∑j=1pAjTAjx−bj,∂f∂Ajx,A=2Ajx−bjxT.Finally, it follows from Section 1.9 in [14] that a norm of ∇fx,A∈ℝn×ℝm×np is defined by
(6)∇fx,A=2∑j=1pAjTAjx−bj22+∑j=1p2Ajx−bjxTfr2,where ∥·∥fr is the Frobenius norm.


## 2.2. The MSEN Problem

Given matricesAj∈ℝm×n and vectors bj∈ℝm, for j=1,⋯,p, the MSEN problem finds a vector x^∈ℝn such that
(7)∑j=1pAjx^−bj2≔minx∑j=1pAjx−bj2.The MSEN problem emerges in many applications, such as the superresolution mapping, VLSL design, transportation and logistics, and the Steiner minimal tree problem. Furthermore, there are several algorithms to estimate optimalx^ in the MSEN problem, for example, the predictor-corrector method developed in [2], the primal-dual algorithm created in [5], and the interior-point technique presented in [6].Note that ifx^ is a solution of problem (7) such that
(8)∑j=1pAjx^−bj2=0,then x^ is a solution of p equivalent linear systems Ajx=bj. However, a solution to the problem (7) does not always satisfy the condition (8) because the accuracy of the MSEN problem is not yet satisfactory. To increase this accuracy, we also need to increase the number of parameters to optimize. Therefore, in the below section, we proposed a new minimization method called the problem of minimizing a sum of squared Euclidean norms with rank constrain, or simply the MSSEN-RC problem.


## 3. A Theoretical Study for the MSSEN-RC Problem

The goal of the MSSEN-RC problem is to estimate a vectorx^ and rank-constrained matrices A^j, for j=1,⋯,p, such that they minimize a sums of squared Euclidean norms. Mathematically, given vectors bj∈ℝm and positive integer rj≤minm,n, for j=1,⋯,p, the MSSEN-RC problem finds x^∈ℝn and A^j∈ℝrjm×n, for j=1,⋯,p, such that
(9)∑j=1pA^jx^−bj22≔minx,A1,⋯,Ap∑j=1pAjx−bj22.In Theorem1 below, we establish that the MSSEN-RC problem (9) has an infinite number of solutions.Theorem 1.
The MSSEN-RC problem has a global minimizer, and it is not unique.Proof.
A solution of the MSSEN-RC problem is given by any nonnull vectorx^∈ℝn and
(10)A^j=bjx∧†+HjIn−x^x∧†,for j=1,⋯,p, where x∧† is the Moore-Penrose inverse of x^ and Hj∈ℝm×m is an arbitrary matrix such that rankA^j≤rj. Based on the fact that x∧†x^=1 (see, e.g., [20]), it is clear that x^ and each A^j given by (10) are a global minimizer of problem (9) because
(11)∑j=1pA^jx^−bj22=∑j=1pbjx∧†+HjIn−x^x∧†x^−bj22=∑j=1pbjx∧†x^+Hjx^−x^x∧†x^−bj22=∑j=1pbj−bj22=0.Finally, solution of problem (9) is not unique because x^ and Hj are arbitrarily chosen.Some relevant remarks on the MSSEN-RC problem are presented below.Remark 2.
Note that ifH is the null matrix in (10), then A^j=bjx∧† is a solution for problem (9), where rankA^j=1, for all j=1,⋯,p. However, finding matrices Hj in (10) such that rankA^j≤rj and rankA~j≠1 is not straightforward. Therefore, in this paper, we consider the MSSEN-RC problem when rj≥2 and 2≤rankA^j≤rj, for all j=1,⋯,p.Remark 3.
Proof of Theorem1 shows the existence of a global solution to problem (9). Additionally, it follows from (11) that if b1,⋯,bp are given, then (8) is true and, therefore, A^jx^=bj, for all j=1,⋯,p.. Thus, we conclude that the MSSEN-RC problem allows us to obtain p equivalent linear systems such that each A^j is rank-constrained.Remark 4.
There are two differences between the MSEN problem and the MSSEN-RC problem. First, the MSSEN-RC problem estimate the optimal vectorx^ and rank-constrained matrices A^j, while the MSEN problem just estimates optimal vector x^. This modification is necessary to obtain p equivalent linear systems. Second, we consider a least-squares problem in (9) instead of a least-absolute problem as in (7). This modification allows us to obtain a stable solution in (9) because the objective function in the MSSEN-RC problem is a differentiable function, unlike the objective function in the MSEN problem.


## 4. A Regularized Alternating Least-Squares Method to Solve (9)

The MSSEN-RC problem can be solved by the alternating least-squares method, i.e., alternatively solving the following two subproblems:(12)xk+1=argminx∈ℝmfx,A1k,⋯,Apk,Ajk+1=argminAj∈ℝrjm×nfxk+1,A1k+1,⋯,Aj−1k+1,Aj,Aj+1k,⋯,Apk,for j=1,⋯,p and k=0,1,2,⋯. Generally, an alternating least-squares method can produce a converging sequence with limit points that are critical points of the problem under certain conditions. Two of these conditions are that each objective function in the alternating minimization method must have a unique solution (see, e.g., Chapter 14 in [14]) or it must be strictly convex (see, e.g., Proposition 5 in [13]). These conditions are not satisfied in (12). Note, for example, that there is no dependency between matrices A1,A2,⋯,Ap in (3) and, therefore, problem (12) can be rewritten as
(13)Ajk+1=argminAj∈ℝrjm×nfjxk+1,Aj=argminAj∈ℝrjm×nAjxk+1−bj22,for j=1,2,⋯,p. A solution of problem (12) is given by (10), i.e.,
(14)Ajk+1=bjxk+1†+HjIn−xk+1xk+1†,where Hj∈ℝm×m is an arbitrary matrix such that rankAjk+1≤rj. Therefore, solution of (12) is not unique.To solve the above issue, we consider regularization terms in (12) that penalize the difference between the previous iterations, which themselves need not be a bounded sequence. Therefore, a regularized alternating least-squares (RALS) method to estimate a solution to the MSSEN-RC problem is given by
(15)xk+1=argminx∈ℝmfx,A1k,⋯,Apk+αkx−xk22,Ajk+1=argminAj∈ℝrjm×nfjxk+1,Aj+βjkAj−Ajkfr2,for j=1,2,⋯,p and k=0,1,2,⋯. Here, αk>0 and βjk>0 are regularization parameters such that αkk=0∞ and βjkk=0∞ are sequences converging to zero. The regularization terms αkx−xk22 and βjkAj−Ajkfr2 are the fitting terms for x and Aj, respectively. The additional regularization terms in (15) penalize the difference between the previous iterates. This approach does not address the degeneracy problem. Moreover, assumptions of convexity and uniqueness of solutions are not required if we are only interested in the critical points (see Section 7 in [13]). The limit points of the RALS algorithm will be shown as the critical points of (3) and not of the regularized version.In Sections4.1 and 4.2 below, we present two algorithms used by the RALS method to estimate a solution to subproblems (15), respectively. The first algorithm is an update to the BFGS method, part of the Quasi-Newton family. The second algorithm uses a generalized low-rank matrix approximation problem to estimate rank-constrained matrices.


### 4.1. A BFGS Method to Solve First Problem in (15)

To solve first problem in (15), we define a function h:ℝn⟶0,∞ such that hx=∑j=1pAjkx−bj22+αkx−xk22, and the minimization problem
(16)minx∈ℝnhx.Note that functionh can be simplified as follows:
hx=∑j=1pAjkx−bj22+αkx−xk22=∑j=1pAjkx−bj22+αkpx−xk22=∑j=1pA~jx−b~j22,whereA~j∈ℝm+n×n and b~j∈ℝm+n such that A~j=Ajkαk/pIn and b~j=bjαk/pxk, for all j=1,2,⋯,p. Therefore, function h is represented by
(17)hx=∑j=1pA~jx−b~j22.Additionally, we consider the gradient function∇h:ℝn⟶ℝn defined by
∇hx=2∑j=1pA~jTA~jx−b~j.Note that ifM=A~1⋮A~p∈ℝmp×nandd=b~1⋮b~p∈ℝmp,then hx=Mx−d22 and, therefore, a solution for the problem (16) is given by the least-squares solution x^=M†d (see, e.g., [20] for more details). Furthermore, we might use the Tikhonov regularization to give preference to a particular solution with desirable properties (see, e.g., [21] for more details). However, the computational complexity increases when m and p are huge in both cases. Moreover, if M is ill-conditioned, then an estimation of M† is not very accurate. Therefore, in this paper, we use an iterative method to compute a solution of problem (16).To estimate a solution to problem (16), we use an update to the classical BFGS method presented by Li and Fukushima in [17]. The classical BFGS method is a well-known Quasi-Newton method for solving unconstrained optimization problems. Moreover, the BFGS method determines the descent direction by preconditioning the gradient with curvature information [22]. Li and Fukushima propose a cautious BFGS update and prove that the objective function to be minimized has Lipschitz continuous gradients. Then, the new BFGS method converges globally using an Armijo-type line search [17]. Note that h in (17) is a Lipschitz continuity function because, for all x,y∈ℝn,
(18)∇hx−∇hy2=2∑j=1pA~jTA~jx−y2≤Lx−y2,where L=2∑j=1pA~jTA~jfr. Last inequality follows from Fact 9.8.41 in [23]. The BFGS method in [24] to estimate a solution of problem (17) is based on the following steps:
(i)
Step 1: choose an initial vectorx¯∈ℝn, initial symmetric and positive-definite matrix D∈ℝn×n and constants ρ,σ∈0,1 and β,ε,tol>0(ii)
Step 2: determine BFGS directionz¯∈ℝn by solving the linear system
Dz¯=−∇hx¯.(iii)
Step 3: find a stepsizeλ that is the largest value in the set ρjj=0,1,⋯ such that the inequality
hx¯+λz¯≤hx¯+σλ∇hx¯Tz¯,is satisfied. This step is the Armijo-type line search.(iv)
Step 4: calculatex¯n=x¯+λz¯(v)
Step 5: computes¯=x¯n−x¯ and y¯=∇hx¯n−∇hx¯(vi)
Step 6: ify¯Ts¯/s¯22≥ε∇hx¯2β, then update D with D−Ds¯s¯TD/s¯TDs¯+y¯y¯T/y¯Ts¯(vii)
Step 7: replacex¯ by x¯n(viii)
Step 8: if∇hx¯2<tol, then quit. Otherwise, go to step 1Algorithm 1: BFGS method.
Input:Ak=A1k,⋯,Apk∈Fp,B=b1⋯bp∈ℝm×p,x¯0∈ℝn,xk∈ℝn,αk>0,tol>0Output:x¯i∈ℝn1 procedure: BFGS_methodAk,B,x¯0,xk,αk,tol2D0=In,ρ=0.2,σ=0.25 and β=∈=0.53  forj=1,2,⋯,p do4Aj~=Ajkαk/pIn and bj~=bjαk/pxk5  fori=0,1,2,⋯ do6    determinez¯i such that Diz¯i=−∇hx¯i7   find smallestj∈0,1,2,⋯ such that λi=ρjhx¯i+λiz¯i≤hx¯i+σλi∇hx¯iTz¯i8xi+1=x¯i+λiz¯i9   if∇hx¯i+12<tol10i=i+111    break12si=x¯i+1−x¯i and y¯i=∇hx¯i+1−∇hx¯i13   ify¯iTs¯i/s¯i22≥ϵ∇hx¯i2β then14Di+1=Di−Dis¯is¯iTDi/s¯iTDis¯i+y¯iy¯iT/y¯iTs¯i15   else16Di+1=Di17  returnx¯iAs mentioned in step 2, initial matrixD is any arbitrary symmetric and positive-definite matrix. Thus, in this paper, we consider as initial matrix the identity matrix of order n, i.e., D=In. Moreover, constants ρ,σ,β,ε can also be chosen arbitrary in the BFGS method. Thus, w.l.g, we use constants ρ=0.2, σ=0.25, and β=ε=0.5. A pseudocode of the BFGS method in [17] is presented in Algorithm 1 below.Some relevant remarks of the BFGS method in Algorithm1 are presented below.Remark 5.
If the level setΩ≔x∈ℝnhx<hx¯0 is bounded, then the BFGS method in [17] converges globally, i.e.,
(19)limi⟶∞inf∇hx¯i2=0.
This result follows from Theorem 3.3 in [17] and the fact that h in (17) is a Lipschitz continuity function (see (18) above). Additionally, Theorem 3.3 in [17] shows that there exists a subsequence of x¯ii=0∞ converging to a stationary point of (16).Remark 6.
In step 7 of Algorithm1, we use the Armijo-type line search to estimate optimal λi. However, Li and Fukushima show in [17] that a Wolfe-type inexact search can also be used in the BFGS method and obtain a global convergence.Remark 7.
Classical BFGS method updates matrixDk+1 using formula given by step 14 of Algorithm 1, even when
(20)y¯iTs¯is¯i22<ε∇hx¯i2β.In this case, if inequality (20) is achieved, then Dk+1 is not necessarily positive-definite, even if Dk is positive-definite. Therefore, Li and Fukushima present in [17] the cautious update shown in steps 13‑16 of Algorithm 1 to ensure the positive-definiteness of Dk+1. This condition is necessary to ensure global convergence. Additionally, the linear system in step 6 of Algorithm 1 can be solved using a fast algorithm for symmetric and positive-definite matrices, for example, the Cholesky decomposition (see, e.g., Theorem 4.18, page 97, in [25]).Remark 8.
In this paper, we choose the BFGS method developed in [17] for its easy computational implementation and high speed and efficiency to estimate a solution of MSEN problem (16). However, there are other methods to estimate a solution of the MSEN problem (see, e.g., [26]).


### 4.2. A Low-Rank Matrix Approximation Problem to Solve Second Problem in (15)

To solve second problem in (15), we consider optimization problem
(21)A^jk+1=argminAj∈ℝrjm×nAjxk+1−bj22+βjkAj−Ajkfr2,for j=1,2,⋯,p. Let Ajxk+1−bj,βjkAj−Ajk∈ℝm×n+1 be a partitioned matrix. Note that
Ajxk+1−bj22+βjkAj−Ajkfr2=Ajxk+1−bj,βjkAj−Ajkfr2=Ajxk+1,βjkAj−bj,βjkAjkfr2=Ajxk+1,βjkIn−bj,βjkAjkfr2.Therefore, ifX¯=xk+1,βjkIn∈ℝn×n+1 and B¯j=bj,βjkAjk∈ℝm×n+1, then problem (21) can be reformulated as follows:
(22)A^jk+1=argminAj∈ℝrjm×nAjX¯−B¯jfr2.A solution of problem (22) was proposed and justified by Friedland and Torokhti in [18], which is given by
(23)A^jk+1=MjrjX¯†,for j=1,2,⋯,p, where Mj=B¯jX¯†X¯ andMjrj is the rj-truncated singular value decomposition (SVD) of Mj, i.e., if Mj=UΣVT is the SVD of Mj, then the rj-truncated SVD of Mj is defined by Mjrj=UrjΣrjVrjT, where Urj and Vrj are formed with the first rj columns of U and V, respectively, and Σrj is formed with the first rj rows and columns of Σ. A pseudocode to solve problem (22), and hence, problem (15), is presented in Algorithm 2 below.Remark 9.
It follows from Theorem 3.1 in [27] that if
(24)rj≤rankX¯B¯jT,then rankA^jk+1=rj, for all j=1,⋯,p. Moreover, it follows from Fact 6.4.2 in [23] that if x and bj are not in the range space of βjk2AjkT and βjk2Ajk, respectively, then
rankX¯B¯jT=rankβjk2AjkT+xk+1bjT=rankβjk2AjkT+1=rj+1.Thus, under the above-mentioned conditions, (24) is true and, therefore, rankA^jk+1=rj≥2, for any rj=2,3,⋯,minm,n. In particular, if we choose 2≤rj≤minm,n, then rankA^jk+1=rj, for all j=1,⋯,p and k=1,2,⋯. Our numerical simulations in Section 5 show that inequality (24) is always true and, therefore, rankA^jk+1=rj≥2, for all j=1,⋯,p and k=0,1,2,⋯.Algorithm 2: Low-rank method.
Input:xk+1∈ℝn, B=b1,⋯,bp∈ℝm\p, Ak=A1k,⋯,Apk∈Fp, r=r_1,⋯,r_p, βk=β1k,⋯,βpkOutput:Ak+1=A1k+1,⋯,Apk+1∈Fp1 procedure: low_rank_methodxk+1,B,Ak,βk,r2 forj=1:p do3X=xk+1,βjkIn4Bj=bj,βjkAjk5Mj=B¯jX¯†X¯6Ajk+1=MjrjX¯†7Ak+1=A1k+1,⋯,Apk+18 returnAk+1


### 4.3. RALS Method and Regularization Parametersαk and βjk

Based on iterative procedure (15) and Algorithms 1 and 2, we show in Algorithm 3 a pseudocode of the RALS method to estimate a solution of the MSSEN-RC problem (9).Algorithm 3: RALS method.
Input:B=b1,⋯,bp∈ℝm×p, r=r_1,⋯,r_px0∈ℝn, A0=A10,⋯,Ap0∈Fp, x¯∈ℝn, tol1>0, tol2>0Output:xk∈ℝn, Ak=A1k,⋯,Apk∈Fp1procedure: RALS_methodB, r, x0, A0, x¯, tol1, tol22k=03  fork=0,1,2,⋯Do4    computeαk5xk+1= BFGS_method Ak, B, x¯, xk, αk, tol16    computeβk=β1k,⋯,βpk7Ak+1=low_rank_method xk+1,B,Ak,βk,r8    if∇fxk+1,Ak+1<tol2, then9k=k+110      break11 returnxk, AkMoreover, recall that the goal in regularization is to produce a solution as close as possible to the accurate answer. Therefore, in this section, we analyze two approaches to compute regularization parametersαk and βjk in steps 4 and 6, respectively, in Algorithm 3.1. Optimal estimation: first, note that optimization problem in (16) can be reformulated as
(25)minz∈ℝnFz−g22+αkz22,where
F=A1k⋮Apk∈ℝmp×n,g=b1−A1kxk⋮bp−Apkxk∈ℝmp,and z=x−xk. A solution for problem (16) is given by x^=z^+xk, where z^ is a solution for problem (25). We consider two cases to compute optimal αk:If the variance of the datag is unknown, we can use the method of generalized cross-validation (GVC) [28] to get an approximate optimal value for αk. The GCV estimate of αk is given by
αk=mpα^,where α^ is the minimizer of
Vα=1/mpKαg221/mptraceKα2,where Kα=Imp−FFTF+mpαImp−1FT.If the variance of the datag is known, then g consists of the true data plus random noise, i.e., g=gt+e, where e∈ℝm is a random vector from a probability distribution with mean 0 and standard deviation s. Let F=UΣVT be the SVD of F, where U=u1⋯um∈ℝm×m and V=v1⋯vn∈ℝn×n are orthogonal matrices, and Σ=diagσ1,⋯,σminm,n∈ℝm×n is a generalized diagonal matrix. It follows from [29] that optimal αk can be estimated by α^, where α^ is the minimizer of
Wα=∑i=1nλi2ασi2+α3−∑i=tnλi2σi2+α2−s2∑i=1t−11σi2+α2,where λi=uiTg and t is a suitable index, depending on the standard deviation s. Other methods to estimate optimal αk, when variance of data is known, are presented in [24, 30].Otherwise, to compute optimalβjk>0 in (21), for j=1,2,⋯,p, we need to solve the problem
minβjk>0A^jk+1−A~jk+1fr2,where A^jk+1 is the solution for regularization problem (21), which is given by (23), and A~jk+1 is a solution for no regularized problem
minAj∈ℝrjm×nAjxk+1−bj22,which is given by (14).2. Approximate estimation: the iterative procedure (15) is an iterative Tikhonov regularization method. As mentioned in [31], this method can be understood in two ways: either αk or βjk is fixed and the iteration number is the main steering parameter or the iteration number is set a priori, and αk and βjk are tuned. If the tuning parameter is based on the parameter choice rule, both procedures give a convergent regularization. The parameter choice must be related to the noise level in the data according to a well-known result by Bakushinskii [32, 33]. Otherwise, a choice independent of the noise level cannot be expected to give convergent results in the worst case. As mentioned in [31], the two schemes above can be mixed to approximate optimal αk and βjk by choosing a different regularization parameter at each iteration, e.g., a geometrically decaying sequence αk=βjk=qk with 0<q<1, for all k=0,1,2,⋯. In fact, this is the strategy that we use in the numerical experiments in Section 5 with q=0.5.


### 4.4. Convergence Analysis

In Theorems10 and 11 below, we present the convergence analysis of proposed RALS method in Algorithm 3. In both theorems, we use the following notations. We define Ak=A1k,⋯,Apk∈Fp and Ajk+1∈Fp such that
Ajk+1=A1k+1,⋯,Aj−1k+1,Ajk+1,Aj+1k,⋯,Apk,for j=1,⋯,p−1, and Apk+1=Ak+1. Let 0n and 0m×n be the zero vector in ℝn and zero matrix in ℝm×n, respectively. We define the zero element 0m×np in ℝm×np as 0m×np=0m×n,⋯,0m×n.Theorem 10.
Letxk,Akk=0∞ be a sequence obtained from Algorithm 3. Then, sequence fxk,Akk=0∞ is convergent.Proof.
Algorithm3 is based on subproblems (15). These subproblems provide the following inequalities:
(26)fxk+1,Ak+αkxk+1−xk22≤fxk,Ak,(27)fxk+1,Ajk+1+βjkAjk+1−Ajkfr2≤fxk+1,Aj−1k+1,for j=1,·⋯,p. It follows from (26) and (27) that
fxk+1,Ak≤fxk,Ak−αkxk+1−xk22≤fxk,Ak,(28)fxk+1,Ajk+1≤fxk+1,Aj−1k+1−βjkAjk+1−Ajkfr2≤fxk+1,Aj−1k+1,and therefore
(29)fxk+1,Ak≤fxk,Ak,fxk+1,Ajk+1≤fxk+1,Aj−1k+1,for all j=1,·⋯,p. Using inequalities (28)-(29) repeatedly, we have
fxk,Ak≥fxk+1,Ak≥fxk+1,A1k+1≥fxk+1,A2k+1⋮≥fxk+1,Ap−1k+1≥fxk+1,Ak+1,and therefore
(30)fxk,Ak≥fxk+1,Ak+1,for all k=0,1,2,⋯. It follows from (30) that fxk,Akk=0∞ is a decreasing sequence. Moreover, fxk,Akk=0∞ is bounded below because fxk,Ak≥0, for all k=0,1,2,⋯. Thus, for the monotone convergence theorem (see, e.g., Theorem 2.3, page 39, in [34]), we conclude that fxk,Akk=0∞ is convergent.Theorem 11.
Letxk,Akk=0∞ be a sequence obtained from Algorithm 3. If xk,Ak⟶x∗,A∗ as k⟶∞, then.
1.fxk,Ak⟶fx∗,A∗ as k⟶∞.
2.x∗,A∗ is a critical point of function fx,A, i.e., ∇fx∗,A∗=0n,0m×np.Proof.
1. From inequality (30), we obtain that
(31)fxk+1,Ak+1≤fxk+1,Ajk+1≤fxk+1,Ak≤fxk,Ak,for all j=1,2,⋯,p. It follows from the fact that xk,Ak⟶x∗,A∗ as k⟶∞, inequalities in (31), the continuity of f and the squeeze theorem that
(32)limk⟶∞fxk+1,Ajk+1=limk⟶∞fxk+1,Ak=limk⟶∞fxk,Ak=fx∗,A∗,for all j=1,⋯,p, i.e., fxk,Ak⟶fx∗,A∗ as k⟶∞.
2. It follows from (32) that making the limits in (26) and (27), we have
limk⟶∞xk+1−xk22=0andlimk⟶∞Ajk+1−Ajkfr2=0,which implies
(33)limk⟶∞xk+1,Ak=x∗,A∗,limk⟶∞xk+1,Ajk+1=x∗,A∗,for j=1,2,⋯,p. From (15), we deduce that xk+1 and Ajk+1 minimize fx,Ak+αkx−xk22 and fxk+1,A1k+1,⋯,Aj−1k+1,Aj,Aj+1k,⋯,Apk+βjkAj−Ajkfr2, respectively. Therefore,
(34)∂f∂xxk+1,Ak+2αkxk+1−xk=0n,∂f∂Ajxk+1,Ajk+1+2βjkAjk+1−Ajk=0m×n,for j=1,2,⋯,p. Then, taking k⟶∞ in (34), using arguments in (33), the continuity of ∇f and the fact that αkk=0∞ and βjkk=0∞ are sequences converging to zero, we obtain
(35)∂f∂xx∗,A∗=0n,∂f∂Ajx∗,A∗=0m×n,for j=1,2,⋯,p. Finally, it follows from (35) that x∗,A∗ is a critical point of function fx,A.
Thus, the proof of Theorem11 is completed.Some remarks of Theorem11 and proposed RALS method in Algorithm 3 are presented below.Remark 12.
As a result of Theorem11, the RALS method solves the same alternating least-squares cost function. Moreover, we proved that the limit point obtained from the RALS algorithm is a critical point of fx,A. However, our numerical simulations in Section 5 show that this critical point is typically also globally optimal.Remark 13.
Theorem11 is a conditional convergence proof, impinging upon the existence of the limit point x∗,A∗. Thus, this result does not address the degeneracy problem. An analysis of the existence of the limit of the RALS method is not included in this work as it is a challenging problem that would require careful study.


## 4.1. A BFGS Method to Solve First Problem in (15)

To solve first problem in (15), we define a function h:ℝn⟶0,∞ such that hx=∑j=1pAjkx−bj22+αkx−xk22, and the minimization problem
(16)minx∈ℝnhx.Note that functionh can be simplified as follows:
hx=∑j=1pAjkx−bj22+αkx−xk22=∑j=1pAjkx−bj22+αkpx−xk22=∑j=1pA~jx−b~j22,whereA~j∈ℝm+n×n and b~j∈ℝm+n such that A~j=Ajkαk/pIn and b~j=bjαk/pxk, for all j=1,2,⋯,p. Therefore, function h is represented by
(17)hx=∑j=1pA~jx−b~j22.Additionally, we consider the gradient function∇h:ℝn⟶ℝn defined by
∇hx=2∑j=1pA~jTA~jx−b~j.Note that ifM=A~1⋮A~p∈ℝmp×nandd=b~1⋮b~p∈ℝmp,then hx=Mx−d22 and, therefore, a solution for the problem (16) is given by the least-squares solution x^=M†d (see, e.g., [20] for more details). Furthermore, we might use the Tikhonov regularization to give preference to a particular solution with desirable properties (see, e.g., [21] for more details). However, the computational complexity increases when m and p are huge in both cases. Moreover, if M is ill-conditioned, then an estimation of M† is not very accurate. Therefore, in this paper, we use an iterative method to compute a solution of problem (16).To estimate a solution to problem (16), we use an update to the classical BFGS method presented by Li and Fukushima in [17]. The classical BFGS method is a well-known Quasi-Newton method for solving unconstrained optimization problems. Moreover, the BFGS method determines the descent direction by preconditioning the gradient with curvature information [22]. Li and Fukushima propose a cautious BFGS update and prove that the objective function to be minimized has Lipschitz continuous gradients. Then, the new BFGS method converges globally using an Armijo-type line search [17]. Note that h in (17) is a Lipschitz continuity function because, for all x,y∈ℝn,
(18)∇hx−∇hy2=2∑j=1pA~jTA~jx−y2≤Lx−y2,where L=2∑j=1pA~jTA~jfr. Last inequality follows from Fact 9.8.41 in [23]. The BFGS method in [24] to estimate a solution of problem (17) is based on the following steps:
(i)
Step 1: choose an initial vectorx¯∈ℝn, initial symmetric and positive-definite matrix D∈ℝn×n and constants ρ,σ∈0,1 and β,ε,tol>0(ii)
Step 2: determine BFGS directionz¯∈ℝn by solving the linear system
Dz¯=−∇hx¯.(iii)
Step 3: find a stepsizeλ that is the largest value in the set ρjj=0,1,⋯ such that the inequality
hx¯+λz¯≤hx¯+σλ∇hx¯Tz¯,is satisfied. This step is the Armijo-type line search.(iv)
Step 4: calculatex¯n=x¯+λz¯(v)
Step 5: computes¯=x¯n−x¯ and y¯=∇hx¯n−∇hx¯(vi)
Step 6: ify¯Ts¯/s¯22≥ε∇hx¯2β, then update D with D−Ds¯s¯TD/s¯TDs¯+y¯y¯T/y¯Ts¯(vii)
Step 7: replacex¯ by x¯n(viii)
Step 8: if∇hx¯2<tol, then quit. Otherwise, go to step 1Algorithm 1: BFGS method.
Input:Ak=A1k,⋯,Apk∈Fp,B=b1⋯bp∈ℝm×p,x¯0∈ℝn,xk∈ℝn,αk>0,tol>0Output:x¯i∈ℝn1 procedure: BFGS_methodAk,B,x¯0,xk,αk,tol2D0=In,ρ=0.2,σ=0.25 and β=∈=0.53  forj=1,2,⋯,p do4Aj~=Ajkαk/pIn and bj~=bjαk/pxk5  fori=0,1,2,⋯ do6    determinez¯i such that Diz¯i=−∇hx¯i7   find smallestj∈0,1,2,⋯ such that λi=ρjhx¯i+λiz¯i≤hx¯i+σλi∇hx¯iTz¯i8xi+1=x¯i+λiz¯i9   if∇hx¯i+12<tol10i=i+111    break12si=x¯i+1−x¯i and y¯i=∇hx¯i+1−∇hx¯i13   ify¯iTs¯i/s¯i22≥ϵ∇hx¯i2β then14Di+1=Di−Dis¯is¯iTDi/s¯iTDis¯i+y¯iy¯iT/y¯iTs¯i15   else16Di+1=Di17  returnx¯iAs mentioned in step 2, initial matrixD is any arbitrary symmetric and positive-definite matrix. Thus, in this paper, we consider as initial matrix the identity matrix of order n, i.e., D=In. Moreover, constants ρ,σ,β,ε can also be chosen arbitrary in the BFGS method. Thus, w.l.g, we use constants ρ=0.2, σ=0.25, and β=ε=0.5. A pseudocode of the BFGS method in [17] is presented in Algorithm 1 below.Some relevant remarks of the BFGS method in Algorithm1 are presented below.Remark 5.
If the level setΩ≔x∈ℝnhx<hx¯0 is bounded, then the BFGS method in [17] converges globally, i.e.,
(19)limi⟶∞inf∇hx¯i2=0.
This result follows from Theorem 3.3 in [17] and the fact that h in (17) is a Lipschitz continuity function (see (18) above). Additionally, Theorem 3.3 in [17] shows that there exists a subsequence of x¯ii=0∞ converging to a stationary point of (16).Remark 6.
In step 7 of Algorithm1, we use the Armijo-type line search to estimate optimal λi. However, Li and Fukushima show in [17] that a Wolfe-type inexact search can also be used in the BFGS method and obtain a global convergence.Remark 7.
Classical BFGS method updates matrixDk+1 using formula given by step 14 of Algorithm 1, even when
(20)y¯iTs¯is¯i22<ε∇hx¯i2β.In this case, if inequality (20) is achieved, then Dk+1 is not necessarily positive-definite, even if Dk is positive-definite. Therefore, Li and Fukushima present in [17] the cautious update shown in steps 13‑16 of Algorithm 1 to ensure the positive-definiteness of Dk+1. This condition is necessary to ensure global convergence. Additionally, the linear system in step 6 of Algorithm 1 can be solved using a fast algorithm for symmetric and positive-definite matrices, for example, the Cholesky decomposition (see, e.g., Theorem 4.18, page 97, in [25]).Remark 8.
In this paper, we choose the BFGS method developed in [17] for its easy computational implementation and high speed and efficiency to estimate a solution of MSEN problem (16). However, there are other methods to estimate a solution of the MSEN problem (see, e.g., [26]).


## 4.2. A Low-Rank Matrix Approximation Problem to Solve Second Problem in (15)

To solve second problem in (15), we consider optimization problem
(21)A^jk+1=argminAj∈ℝrjm×nAjxk+1−bj22+βjkAj−Ajkfr2,for j=1,2,⋯,p. Let Ajxk+1−bj,βjkAj−Ajk∈ℝm×n+1 be a partitioned matrix. Note that
Ajxk+1−bj22+βjkAj−Ajkfr2=Ajxk+1−bj,βjkAj−Ajkfr2=Ajxk+1,βjkAj−bj,βjkAjkfr2=Ajxk+1,βjkIn−bj,βjkAjkfr2.Therefore, ifX¯=xk+1,βjkIn∈ℝn×n+1 and B¯j=bj,βjkAjk∈ℝm×n+1, then problem (21) can be reformulated as follows:
(22)A^jk+1=argminAj∈ℝrjm×nAjX¯−B¯jfr2.A solution of problem (22) was proposed and justified by Friedland and Torokhti in [18], which is given by
(23)A^jk+1=MjrjX¯†,for j=1,2,⋯,p, where Mj=B¯jX¯†X¯ andMjrj is the rj-truncated singular value decomposition (SVD) of Mj, i.e., if Mj=UΣVT is the SVD of Mj, then the rj-truncated SVD of Mj is defined by Mjrj=UrjΣrjVrjT, where Urj and Vrj are formed with the first rj columns of U and V, respectively, and Σrj is formed with the first rj rows and columns of Σ. A pseudocode to solve problem (22), and hence, problem (15), is presented in Algorithm 2 below.Remark 9.
It follows from Theorem 3.1 in [27] that if
(24)rj≤rankX¯B¯jT,then rankA^jk+1=rj, for all j=1,⋯,p. Moreover, it follows from Fact 6.4.2 in [23] that if x and bj are not in the range space of βjk2AjkT and βjk2Ajk, respectively, then
rankX¯B¯jT=rankβjk2AjkT+xk+1bjT=rankβjk2AjkT+1=rj+1.Thus, under the above-mentioned conditions, (24) is true and, therefore, rankA^jk+1=rj≥2, for any rj=2,3,⋯,minm,n. In particular, if we choose 2≤rj≤minm,n, then rankA^jk+1=rj, for all j=1,⋯,p and k=1,2,⋯. Our numerical simulations in Section 5 show that inequality (24) is always true and, therefore, rankA^jk+1=rj≥2, for all j=1,⋯,p and k=0,1,2,⋯.Algorithm 2: Low-rank method.
Input:xk+1∈ℝn, B=b1,⋯,bp∈ℝm\p, Ak=A1k,⋯,Apk∈Fp, r=r_1,⋯,r_p, βk=β1k,⋯,βpkOutput:Ak+1=A1k+1,⋯,Apk+1∈Fp1 procedure: low_rank_methodxk+1,B,Ak,βk,r2 forj=1:p do3X=xk+1,βjkIn4Bj=bj,βjkAjk5Mj=B¯jX¯†X¯6Ajk+1=MjrjX¯†7Ak+1=A1k+1,⋯,Apk+18 returnAk+1


## 4.3. RALS Method and Regularization Parametersαk and βjk

Based on iterative procedure (15) and Algorithms 1 and 2, we show in Algorithm 3 a pseudocode of the RALS method to estimate a solution of the MSSEN-RC problem (9).Algorithm 3: RALS method.
Input:B=b1,⋯,bp∈ℝm×p, r=r_1,⋯,r_px0∈ℝn, A0=A10,⋯,Ap0∈Fp, x¯∈ℝn, tol1>0, tol2>0Output:xk∈ℝn, Ak=A1k,⋯,Apk∈Fp1procedure: RALS_methodB, r, x0, A0, x¯, tol1, tol22k=03  fork=0,1,2,⋯Do4    computeαk5xk+1= BFGS_method Ak, B, x¯, xk, αk, tol16    computeβk=β1k,⋯,βpk7Ak+1=low_rank_method xk+1,B,Ak,βk,r8    if∇fxk+1,Ak+1<tol2, then9k=k+110      break11 returnxk, AkMoreover, recall that the goal in regularization is to produce a solution as close as possible to the accurate answer. Therefore, in this section, we analyze two approaches to compute regularization parametersαk and βjk in steps 4 and 6, respectively, in Algorithm 3.1. Optimal estimation: first, note that optimization problem in (16) can be reformulated as
(25)minz∈ℝnFz−g22+αkz22,where
F=A1k⋮Apk∈ℝmp×n,g=b1−A1kxk⋮bp−Apkxk∈ℝmp,and z=x−xk. A solution for problem (16) is given by x^=z^+xk, where z^ is a solution for problem (25). We consider two cases to compute optimal αk:If the variance of the datag is unknown, we can use the method of generalized cross-validation (GVC) [28] to get an approximate optimal value for αk. The GCV estimate of αk is given by
αk=mpα^,where α^ is the minimizer of
Vα=1/mpKαg221/mptraceKα2,where Kα=Imp−FFTF+mpαImp−1FT.If the variance of the datag is known, then g consists of the true data plus random noise, i.e., g=gt+e, where e∈ℝm is a random vector from a probability distribution with mean 0 and standard deviation s. Let F=UΣVT be the SVD of F, where U=u1⋯um∈ℝm×m and V=v1⋯vn∈ℝn×n are orthogonal matrices, and Σ=diagσ1,⋯,σminm,n∈ℝm×n is a generalized diagonal matrix. It follows from [29] that optimal αk can be estimated by α^, where α^ is the minimizer of
Wα=∑i=1nλi2ασi2+α3−∑i=tnλi2σi2+α2−s2∑i=1t−11σi2+α2,where λi=uiTg and t is a suitable index, depending on the standard deviation s. Other methods to estimate optimal αk, when variance of data is known, are presented in [24, 30].Otherwise, to compute optimalβjk>0 in (21), for j=1,2,⋯,p, we need to solve the problem
minβjk>0A^jk+1−A~jk+1fr2,where A^jk+1 is the solution for regularization problem (21), which is given by (23), and A~jk+1 is a solution for no regularized problem
minAj∈ℝrjm×nAjxk+1−bj22,which is given by (14).2. Approximate estimation: the iterative procedure (15) is an iterative Tikhonov regularization method. As mentioned in [31], this method can be understood in two ways: either αk or βjk is fixed and the iteration number is the main steering parameter or the iteration number is set a priori, and αk and βjk are tuned. If the tuning parameter is based on the parameter choice rule, both procedures give a convergent regularization. The parameter choice must be related to the noise level in the data according to a well-known result by Bakushinskii [32, 33]. Otherwise, a choice independent of the noise level cannot be expected to give convergent results in the worst case. As mentioned in [31], the two schemes above can be mixed to approximate optimal αk and βjk by choosing a different regularization parameter at each iteration, e.g., a geometrically decaying sequence αk=βjk=qk with 0<q<1, for all k=0,1,2,⋯. In fact, this is the strategy that we use in the numerical experiments in Section 5 with q=0.5.


## 4.4. Convergence Analysis

In Theorems10 and 11 below, we present the convergence analysis of proposed RALS method in Algorithm 3. In both theorems, we use the following notations. We define Ak=A1k,⋯,Apk∈Fp and Ajk+1∈Fp such that
Ajk+1=A1k+1,⋯,Aj−1k+1,Ajk+1,Aj+1k,⋯,Apk,for j=1,⋯,p−1, and Apk+1=Ak+1. Let 0n and 0m×n be the zero vector in ℝn and zero matrix in ℝm×n, respectively. We define the zero element 0m×np in ℝm×np as 0m×np=0m×n,⋯,0m×n.Theorem 10.
Letxk,Akk=0∞ be a sequence obtained from Algorithm 3. Then, sequence fxk,Akk=0∞ is convergent.Proof.
Algorithm3 is based on subproblems (15). These subproblems provide the following inequalities:
(26)fxk+1,Ak+αkxk+1−xk22≤fxk,Ak,(27)fxk+1,Ajk+1+βjkAjk+1−Ajkfr2≤fxk+1,Aj−1k+1,for j=1,·⋯,p. It follows from (26) and (27) that
fxk+1,Ak≤fxk,Ak−αkxk+1−xk22≤fxk,Ak,(28)fxk+1,Ajk+1≤fxk+1,Aj−1k+1−βjkAjk+1−Ajkfr2≤fxk+1,Aj−1k+1,and therefore
(29)fxk+1,Ak≤fxk,Ak,fxk+1,Ajk+1≤fxk+1,Aj−1k+1,for all j=1,·⋯,p. Using inequalities (28)-(29) repeatedly, we have
fxk,Ak≥fxk+1,Ak≥fxk+1,A1k+1≥fxk+1,A2k+1⋮≥fxk+1,Ap−1k+1≥fxk+1,Ak+1,and therefore
(30)fxk,Ak≥fxk+1,Ak+1,for all k=0,1,2,⋯. It follows from (30) that fxk,Akk=0∞ is a decreasing sequence. Moreover, fxk,Akk=0∞ is bounded below because fxk,Ak≥0, for all k=0,1,2,⋯. Thus, for the monotone convergence theorem (see, e.g., Theorem 2.3, page 39, in [34]), we conclude that fxk,Akk=0∞ is convergent.Theorem 11.
Letxk,Akk=0∞ be a sequence obtained from Algorithm 3. If xk,Ak⟶x∗,A∗ as k⟶∞, then.
1.fxk,Ak⟶fx∗,A∗ as k⟶∞.
2.x∗,A∗ is a critical point of function fx,A, i.e., ∇fx∗,A∗=0n,0m×np.Proof.
1. From inequality (30), we obtain that
(31)fxk+1,Ak+1≤fxk+1,Ajk+1≤fxk+1,Ak≤fxk,Ak,for all j=1,2,⋯,p. It follows from the fact that xk,Ak⟶x∗,A∗ as k⟶∞, inequalities in (31), the continuity of f and the squeeze theorem that
(32)limk⟶∞fxk+1,Ajk+1=limk⟶∞fxk+1,Ak=limk⟶∞fxk,Ak=fx∗,A∗,for all j=1,⋯,p, i.e., fxk,Ak⟶fx∗,A∗ as k⟶∞.
2. It follows from (32) that making the limits in (26) and (27), we have
limk⟶∞xk+1−xk22=0andlimk⟶∞Ajk+1−Ajkfr2=0,which implies
(33)limk⟶∞xk+1,Ak=x∗,A∗,limk⟶∞xk+1,Ajk+1=x∗,A∗,for j=1,2,⋯,p. From (15), we deduce that xk+1 and Ajk+1 minimize fx,Ak+αkx−xk22 and fxk+1,A1k+1,⋯,Aj−1k+1,Aj,Aj+1k,⋯,Apk+βjkAj−Ajkfr2, respectively. Therefore,
(34)∂f∂xxk+1,Ak+2αkxk+1−xk=0n,∂f∂Ajxk+1,Ajk+1+2βjkAjk+1−Ajk=0m×n,for j=1,2,⋯,p. Then, taking k⟶∞ in (34), using arguments in (33), the continuity of ∇f and the fact that αkk=0∞ and βjkk=0∞ are sequences converging to zero, we obtain
(35)∂f∂xx∗,A∗=0n,∂f∂Ajx∗,A∗=0m×n,for j=1,2,⋯,p. Finally, it follows from (35) that x∗,A∗ is a critical point of function fx,A.
Thus, the proof of Theorem11 is completed.Some remarks of Theorem11 and proposed RALS method in Algorithm 3 are presented below.Remark 12.
As a result of Theorem11, the RALS method solves the same alternating least-squares cost function. Moreover, we proved that the limit point obtained from the RALS algorithm is a critical point of fx,A. However, our numerical simulations in Section 5 show that this critical point is typically also globally optimal.Remark 13.
Theorem11 is a conditional convergence proof, impinging upon the existence of the limit point x∗,A∗. Thus, this result does not address the degeneracy problem. An analysis of the existence of the limit of the RALS method is not included in this work as it is a challenging problem that would require careful study.


## 5. Numerical Experiments

In this section, we show some numerical experiments of the RALS method in Algorithm3 to estimate a solution to the MSSEN-RC problem (9). Algorithm 3 was implemented in GNU Octave 6.1.0. The numerical experiments were run on a desktop computer with 2.20 GHz processor (Intel Xeon E5 2660 v2) and 48.00 GB RAM.As mentioned in Section1, the proposed MSSEN-RC problem in (9) is a new optimization problem. For that reason, we cannot compare Algorithm 3 with another method that solves that problem. However, the MSSEN-RC problem generalizes the classical MSEN problem given in (7). In the following experiments, we compare Algorithm 3 with a predictor-corrector method given in [2], which estimates a solution of the MSEN problem. This predictor-corrector technique is derived from a primal-dual interior-point algorithm by applying Newton’s method directly to a system of nonlinear equations characterizing primal and dual feasibility and a perturbed complementary condition. The stopping condition of this method is given by the duality gap, i.e.,
gapx=∑i=1pzj2−yjTzj,where yj∈ℝm is a dual variable from the MSEN problem and zj=bj−Ajx, for all j=1,2,⋯,p (more details of the predictor-corrector algorithm are given in [2]). Based on Algorithm 3 and the predictor-corrector method, we numerically show that the minimum obtained in the MSSEN-RC problem is less than that of the MSEN problem, i.e.,
minx,Afx,A≤minxgx,where fx,A=fx,A1,⋯,Ap=∑j=1pAjx−bj22 and gx=∑j=1pAjx−bj2.


### 5.1. Experiment 1

We consider the first six numerical examples given in Section 7 in [35]. These examples are used to solve the MSEN problem (9). Therefore, we generalized the numerical results in [35] to the proposed MSSEN-RC problem (9) and Algorithm 3. We defined x0∈ℝ2, bj∈ℝ2, and Aj0∈ℝ2×2, for j=1,2,3, such that b1=−1,0T, b2=0,σT, b3=1,0T, A10=A30=I2 and A20=σI2, where σ>0. In Algorithm 3, we use r1=r2=r3=2, tol1=tol2=10−10, x¯=1,1T, and αk=βjk=0.5k.Our numerical results are summarized in Tables1, 2, 3, and 4. In Tables 1 and 2, k is the number of iterations. In Tables 3 and 4, we present the solutions to each numerical simulation given by Algorithm 3 and the predictor-corrector method in [2], respectively. Additionally, Figure 1 shows diagrams of the errors associated with the numerical simulations in Table 1.Table 1
Experiment 1: results of problem (9), given by Algorithm 3.
Simulationx0σk∇fxk,Akfxk,Aka3,2T265.6034×10−121.1874×10−23b1,10−6T267.0912×10−121.9943×10−23c1.000001,−10−6T267.0911×10−121.9943×10−23d1.001,−10−3T267.0912×10−121.9943×10−23e3,2T163.0609×10−114.8636×10−22f3,2T261.4519×10−119.6208×10−23Table 2
Experiment 1: results of problem (7), given by Algorithm in [2].
Simulationx0σkgapxkgxka3,2T2143.2512×10−114b1,10−6T2246.5823×10−114c1.000001,−10−6T2229.8461×10−114d1.001,−10−3T2217.2201×10−114e3,2T1122.6121×10−112.8453f3,2T2192.9711×10−114Table 3
Experiment 1: solution of problem (9), given by Algorithm 3.
Simulationx∗A1∗A2∗A3∗a0.0390.8080.92−1.28−0.060.0021.99−0.090.032.471.071.18−0.060.003b,c,d0.0130.7910.97−1.28−0.0201.99−0.030.012.531.021.25−0.020e0.1220.6830.68−1.59−0.190.030.97−0.170.101.451.251.24−0.190.03f0.0690.7360.84−1.43−0.1040.011.400.130.061.911.131.25−0.100.01Table 4
Experiment 1: solution of problem (7), given by Algorithm in [2].
Simulationx∗a,b,c,d,f01e0√3/3Figure 1
Diagrams of error∇fxk,Ak versus iterations k of experiment 5.1.
(a)(b)(c)(d)(e)(f)The results listed in Tables1 and 2 show that Algorithm 3 is extremely promising. Moreover, it follows from Table 1 and Figure 1 that Algorithm 3 is able to approximate numerical solutions for all the experiments within reasonable iterations. Particularly, Figure 1 illustrates that the error decreases to zero as the iterations increase, i.e., ∇fxk,Ak⟶0 as k⟶∞. The above sentence follows from Theorem 11.Futhermore, note that matricesA1∗, A2∗, and A3∗ in Table 3 are full-rank, i.e., rankA1∗=rankA2∗=rankA3∗=2. This result follows from inequality (19) in Remark 9, which is true in this experiment, for all j=1,2,3 and k=0,1,2,3,4,5.Finally, it follows from Tables1 and 2 that fxk,Ak<gxk, i.e., the minimum obtained in the MSSEN-RC problem is less than that of the MSEN problem.


### 5.2. Experiment 2

We evaluate the efficiency of Algorithm3 to approximate a solution of the MSSEN-RC problem (9), using random matrices. Specifically, for j=1,⋯,p, we consider matrices Aj0∈ℝrjm×n and vectors bj∈ℝm generated from a normal distribution with zero mean and standard deviation 1. In Algorithm 3, we use tol1=tol2=10−10, x¯∈ℝn is the vector of ones, rj=roundminm,n/2, for j=1,⋯,p and αk=βjk=0.5k.Tables5 and 6 summarizes the numerical results obtained. In both tables, m, n, and p specify the problem dimensions, k denotes the number of iterations, and rj is the rank constraint of each Aj, for j=1,⋯,p. Figure 2 shows diagrams of the errors associated with the numerical simulation in Table 5. Similarly to experiment 5.1, the results presented in Table 5 show that Algorithm 3 is promising. Moreover, Figure 2 shows that the error decreases to zero as the iterations increase, i.e., ∥∇fxk,Ak∥⟶0 as k⟶∞. The above sentence follows from Theorem 11. Based on the numerical data in Table 5 and Figure 2, Algorithm 3 approximates a solution for all the experiments using few iterations.Table 5
Experiment 2: results of problem (9), given by Algorithm 3.
Simulationmnrjpk∇fxk,Akfxk,Aka10205571.0206×10−131.8283×10−27b301581077.9621×10−144.5197×10−28c2040106073.4293×10−122.8296×10−24d406020100103.7677×10−129.3883×10−24Table 6
Experiment 2: results of problem (7), given by Algorithm in [2].
Simulationmnpkgapxkgxka10205201.99×10−1130.55b301510193.29×10−112.55×102c204060211.87×10−111.16×103d4060100171.81×10−114.03×103Figure 2
Diagrams of error∇fxk,Ak versus iterations k of experiment 5.2.
(a)(b)(c)(d)Moreover, we obtain thatrankAjk=rj, for each numerical simulation in this experiment. This result follows from inequality (24) in Remark 9, which is true for all j=1,⋯p.Finally, it follows from Tables5 and 6 that fxk,Ak<gxk, i.e., the minimum obtained in the MSSEN-RC problem is less than that of the MSEN problem.


## 5.1. Experiment 1

We consider the first six numerical examples given in Section 7 in [35]. These examples are used to solve the MSEN problem (9). Therefore, we generalized the numerical results in [35] to the proposed MSSEN-RC problem (9) and Algorithm 3. We defined x0∈ℝ2, bj∈ℝ2, and Aj0∈ℝ2×2, for j=1,2,3, such that b1=−1,0T, b2=0,σT, b3=1,0T, A10=A30=I2 and A20=σI2, where σ>0. In Algorithm 3, we use r1=r2=r3=2, tol1=tol2=10−10, x¯=1,1T, and αk=βjk=0.5k.Our numerical results are summarized in Tables1, 2, 3, and 4. In Tables 1 and 2, k is the number of iterations. In Tables 3 and 4, we present the solutions to each numerical simulation given by Algorithm 3 and the predictor-corrector method in [2], respectively. Additionally, Figure 1 shows diagrams of the errors associated with the numerical simulations in Table 1.Table 1
Experiment 1: results of problem (9), given by Algorithm 3.
Simulationx0σk∇fxk,Akfxk,Aka3,2T265.6034×10−121.1874×10−23b1,10−6T267.0912×10−121.9943×10−23c1.000001,−10−6T267.0911×10−121.9943×10−23d1.001,−10−3T267.0912×10−121.9943×10−23e3,2T163.0609×10−114.8636×10−22f3,2T261.4519×10−119.6208×10−23Table 2
Experiment 1: results of problem (7), given by Algorithm in [2].
Simulationx0σkgapxkgxka3,2T2143.2512×10−114b1,10−6T2246.5823×10−114c1.000001,−10−6T2229.8461×10−114d1.001,−10−3T2217.2201×10−114e3,2T1122.6121×10−112.8453f3,2T2192.9711×10−114Table 3
Experiment 1: solution of problem (9), given by Algorithm 3.
Simulationx∗A1∗A2∗A3∗a0.0390.8080.92−1.28−0.060.0021.99−0.090.032.471.071.18−0.060.003b,c,d0.0130.7910.97−1.28−0.0201.99−0.030.012.531.021.25−0.020e0.1220.6830.68−1.59−0.190.030.97−0.170.101.451.251.24−0.190.03f0.0690.7360.84−1.43−0.1040.011.400.130.061.911.131.25−0.100.01Table 4
Experiment 1: solution of problem (7), given by Algorithm in [2].
Simulationx∗a,b,c,d,f01e0√3/3Figure 1
Diagrams of error∇fxk,Ak versus iterations k of experiment 5.1.
(a)(b)(c)(d)(e)(f)The results listed in Tables1 and 2 show that Algorithm 3 is extremely promising. Moreover, it follows from Table 1 and Figure 1 that Algorithm 3 is able to approximate numerical solutions for all the experiments within reasonable iterations. Particularly, Figure 1 illustrates that the error decreases to zero as the iterations increase, i.e., ∇fxk,Ak⟶0 as k⟶∞. The above sentence follows from Theorem 11.Futhermore, note that matricesA1∗, A2∗, and A3∗ in Table 3 are full-rank, i.e., rankA1∗=rankA2∗=rankA3∗=2. This result follows from inequality (19) in Remark 9, which is true in this experiment, for all j=1,2,3 and k=0,1,2,3,4,5.Finally, it follows from Tables1 and 2 that fxk,Ak<gxk, i.e., the minimum obtained in the MSSEN-RC problem is less than that of the MSEN problem.


## 5.2. Experiment 2

We evaluate the efficiency of Algorithm3 to approximate a solution of the MSSEN-RC problem (9), using random matrices. Specifically, for j=1,⋯,p, we consider matrices Aj0∈ℝrjm×n and vectors bj∈ℝm generated from a normal distribution with zero mean and standard deviation 1. In Algorithm 3, we use tol1=tol2=10−10, x¯∈ℝn is the vector of ones, rj=roundminm,n/2, for j=1,⋯,p and αk=βjk=0.5k.Tables5 and 6 summarizes the numerical results obtained. In both tables, m, n, and p specify the problem dimensions, k denotes the number of iterations, and rj is the rank constraint of each Aj, for j=1,⋯,p. Figure 2 shows diagrams of the errors associated with the numerical simulation in Table 5. Similarly to experiment 5.1, the results presented in Table 5 show that Algorithm 3 is promising. Moreover, Figure 2 shows that the error decreases to zero as the iterations increase, i.e., ∥∇fxk,Ak∥⟶0 as k⟶∞. The above sentence follows from Theorem 11. Based on the numerical data in Table 5 and Figure 2, Algorithm 3 approximates a solution for all the experiments using few iterations.Table 5
Experiment 2: results of problem (9), given by Algorithm 3.
Simulationmnrjpk∇fxk,Akfxk,Aka10205571.0206×10−131.8283×10−27b301581077.9621×10−144.5197×10−28c2040106073.4293×10−122.8296×10−24d406020100103.7677×10−129.3883×10−24Table 6
Experiment 2: results of problem (7), given by Algorithm in [2].
Simulationmnpkgapxkgxka10205201.99×10−1130.55b301510193.29×10−112.55×102c204060211.87×10−111.16×103d4060100171.81×10−114.03×103Figure 2
Diagrams of error∇fxk,Ak versus iterations k of experiment 5.2.
(a)(b)(c)(d)Moreover, we obtain thatrankAjk=rj, for each numerical simulation in this experiment. This result follows from inequality (24) in Remark 9, which is true for all j=1,⋯p.Finally, it follows from Tables5 and 6 that fxk,Ak<gxk, i.e., the minimum obtained in the MSSEN-RC problem is less than that of the MSEN problem.


## 6. Conclusion

In this work, we proposed modifying the classical MSEN problem, which we call the MSSEN-RC problem. The objective of the MSSEN-RC problem is to obtainx^∈ℝn and A^j∈ℝrjm×n, for j=1,⋯,p, such that they minimize function f in (3). The proposed MSSEN-RC problem allows the creation of p equivalent linear systems Ajx=bj, such that each Aj is a rank-constrained matrix. Additionally, we developed an iterative method in Algorithm 3 to estimate a solution to the MSSEN-RC problem. Algorithm 3 is based on a regularized alternating least-squares (RALS) method. We considered two numerical techniques in the RALS method: the BFGS method, part of the Quasi-Newton family, and a generalized low-rank matrix approximation problem to estimate rank-constrained matrices. In Theorem 11, we showed that, under some conditions, Algorithm 3 converges to a critical point of function f. The theoretical development and numerical experiments in this paper demonstrate the advantages of the proposed RALS method in Algorithm 3. As for future work, we will investigate applying the MSSEN-RC problem to solve a real-world inverse problem.



---

*Source: 4838182-2022-05-17.xml*