Iterative methods for solving factorized linear systems
Anna Ma Deanna Needell Aaditya Ramdas


Control Allocation for Wide Area Coordinated Damping
    M. Ehsan Raoufat, Student Member, IEEE, 
		Kevin Tomsovic, Fellow, IEEE, and Seddik M.  Djouadi, Member, IEEE
		
		This work was supported in part by the National Science Foundation under grant No CNS-1239366, and in part by the Engineering Research Center Program of the National Science Foundation and the Department of Energy under NSF Award Number EEC-1041877 and the CURENT Industry Partnership Program.
		
		M. Ehsan Raoufat, Kevin Tomsovic and Seddik M. Djouadi are with the Min H. Kao Department of Electrical Engineering and Computer Science, The University of Tennessee, Knoxville, TN 37996 USA (e-mail: mraoufat@utk.edu).
    A. Mahmoodzadeha.mahmoodzadeh@iau-boukan.ac.ir, B. MalekolkalamiB.Malakolkalami@uok.ac.ir 
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




Stochastic iterative algorithms such as the Kaczmarz and Gauss-Seidel methods have gained recent attention because of their speed, simplicity, and the ability to approximately solve large-scale linear systems of equations without needing to access the entire matrix.  In this work, we consider the setting where we wish to solve a linear system in a large matrix  that is stored in a factorized form, =; this setting either arises naturally in many applications or may be imposed when working with large low-rank datasets for reasons of space required for storage. We propose a variant of the randomized Kaczmarz method for such systems that takes advantage of the factored form, and avoids computing .  We prove an exponential convergence rate and supplement our theoretical guarantees with experimental evidence demonstrating that the factored variant yields significant acceleration in convergence.








§ INTRODUCTION


  Recently, revived interest in stochastic iterative methods like the Kaczmarz <cit.> and Gauss-Seidel <cit.> methods has grown due to the need for large-scale approaches for solving linear systems of equations.  Such methods utilize simple projections and require access to only a single row in a given iteration, hence having a low memory footprint.  For this reason, they are very efficient and practical for solving extremely large, usually highly overdetermined, linear systems.  In this work, we consider algorithms for solving linear systems when the matrix is available in a factorized form.  As we discuss below, such a factorization may arise naturally in the application, or may be constructed explicitly for efficient storage and computation.  We seek a solution to the original system directly from its factorized form, without the need to perform matrix multiplication.

	To that end, borrowing the notation of linear regression from statistics, suppose we want to solve the linear system = with ∈ℂ^m × n.
	However, instead of the full system , we only have access to , such that =. 
	In this case, we want to solve the linear system: 
		
    = ,

	where ∈ℂ^m × k and ∈ℂ^k × n. 
	Instead of taking the product of  and , to form , which may not be desirable, we approach this problem using stochastic iterative methods to solve the individual subsystems
		
    = 
       = ,

in an alternating fashion.
	Note that  in (<ref>) is the vector of unknowns that we want to solve for in (<ref>) and  in (<ref>) is the known right hand side vector of (<ref>). If we substitute (<ref>) into (<ref>), we acquire the full linear system (<ref>). We will often refer to (<ref>) as the “full system” and (<ref>) and (<ref>) as “subsystems”, and say that a system is consistent if it has at least one solution (and inconsistent otherwise). 
	
There are some situations when approximately knowing  would suffice. We assume that (for reasons of interpretability, or for downstream usage) the scientist is genuinely interested in solving the full system, i.e. she is interested in the vector , not in . 

It is arguably of practical interest to give special importance to the case of k < min(m,n), which arises in modern data science as motivated by  the following examples, but we discuss other settings later.

	

 §.§ Motivation
 
		If  is large and low-rank, one may have many reasons to work with a factorization of . We shall discuss three reasons below — algorithmic, infrastructural, and statistical.
		
		Consider data matrices encountered in “recommender systems” in machine learning <cit.>. For concreteness, consider the Netflix (or Amazon, or Yelp) problem, where one has a users-by-movies matrix whose entries correspond to ratings given by users  to movies.  is usually quite well approximated by low-rank matrices — intuitively, many rows and columns are redundant because every row is usually similar to many other rows (corresponding to users with similar tastes), and every column is usually similar to many other columns (corresponding to similar quality movies in the same genre). Usually we have observed only a few entries of , and wish to infer the unseen ratings in order to provide recommendations to different users based on their tastes. Algorithms for “low-rank matrix completion” have proved  to be quite successful in the applied and theoretical machine learning community <cit.>. One popular algorithm,  alternating-minimization <cit.>, chooses a (small) target rank k, and tries to find , such that _ij≈ ()_ij for all the observed entries (i,j) of . As its name suggests, the algorithm alternates between solving for  keeping  fixed and then solving for  keeping  fixed. In this case, at no point does the algorithm even form the entire completed (inferred) matrix , and the algorithm only has access to factors , simply due to algorithmic choices. 

There may be other instances where a data scientist may have access to the full matrix , but in order to reduce the memory storage footprint, or to communicate the data, may explicitly choose to decompose ≈ and discard  to work with the smaller matrices instead. 

		Consider an example motivated by “topic modeling” of text data. 
		Suppose Google has scraped the internet for English documents (or maybe a subset of documents like news articles), to form a document-by-word matrix , where each entry of the matrix indicates the number of times that word occurred in that document. Since many documents could be quite similar in their content (like articles about the same incident covered by different newspapers), this matrix is easily seen to be low-rank. This is a classic setting for applying a machine learning technique called “non-negative matrix factorization” <cit.>, where one decomposes  as the product of two low-rank non-negative matrices ,; the non-negativity is imposed for human interpretability, so that  can be interpreted as a documents-by-topics matrix, and  as a topics-by-words. In this case, we do not have access to  as a result of systems infrastructure constraints (memory/storage/communication).

Often, even for modestly sized data matrices, the relevant “signal” is contained in the leading singular vectors corresponding to large singular values, and the tail of small singular values is often deemed to be “noise”. This is precisely the idea behind the classical topic of principal component analysis (PCA), and the modern machine learning literature has proposed and analyzed a variety of algorithms to approximate the top k left and right singular vectors in a streaming/stochastic/online fashion <cit.>. Hence, the factorization may arise from a purely statistical motivation.

Given a vector  (representing age, or document popularity, for example), suppose the data scientist is interested in regressing  onto , for the purpose of scientific understanding or to take future actions. Can we utilize the available factorization efficiently, designing methods that work directly on the lower dimensional factors  and  rather than computing the full system ?

		
		Our goal will be to propose iterative methods that work directly on the factored system, eliminating the need for a full matrix product and potentially saving computations on the much larger full system.
	
	

 §.§ Main contribution

		We propose two stochastic iterative methods for solving system (<ref>) without computing the product of  and . 
		Both methods utilize iterates of well studied algorithms for solving linear systems.
		When the full system is consistent, the first method, called RK-RK, interlaces iterates of the Randomized Kaczmarz (RK) algorithm to solve each subsystem and finds the optimal solution. 
		When the full system is inconsistent, we introduce the REK-RK method, an interlacing of Randomized Extended Kaczmarz (REK) iterates to solve (<ref>) and RK iterates to solve (<ref>), that converges to the so-called ordinary least squares solution.  We prove linear (“exponential”) convergence to the solution in both cases.
		
	

 §.§ Outline

		In the next section, we provide background and discuss existing work on stochastic methods that solve linear systems. 
		In particular, we describe the RK and REK algorithms as well as the Randomized Gauss-Seidel (RGS) and Randomized Extended Gauss-Seidel (REGS) algorithms. 
		In Section <ref> we investigate variations of settings for subsystems (<ref>) and (<ref>) that arise depending on the consistency and size of .
		Section <ref> introduces our proposed methods, RK-RK and REK-RK. We provide theory that shows linear convergence in expectation to the optimal solution for both methods.
		Finally, we present experiments in Section <ref> and conclude with final remarks and future work in Section <ref>.
		
		

 §.§ Notation

		Here and throughout the paper, matrices and vectors are denoted with boldface letters (uppercase for matrices and lowercase for vectors). 
		We call  the i^th row of the matrix  and  the j^th column of .
		The euclidean norm is denoted by ·_2 and the Frobenius norm by  ·_F.  Lastly, ^* denotes the adjoint (conjugate transpose) of the matrix . Motivated by applications, we allow  to be rank deficient and assume that  and  are full rank.






§ BACKGROUND AND EXISTING WORK

		In this section we summarize existing work on stochastic iterative methods and different variations of linear systems. 
		
		

 §.§ Linear Systems

		Linear systems take on one of three settings determined by the size of the system, rank of the matrix , and the existence of a solution. First we discuss solutions to systems with full rank matrices  then remark on how rank deficiency affects the desired solution.

		 In the full rank underdetermined case, m < n and the system has infinitely many solutions; here, we often want to find the least Euclidean norm solution to (<ref>):
		
    := ^*(^*)^-1 .

		Clearly, =, and all other solutions to an underdetermined system can be written as =  + where = 0.
		
		In the overdetermined setting, we have m > n and the system can have a unique (exact) solution or no solution. If there is a unique solution, the linear system is called an overdetermined consistent system.
		When  is full rank, the optimal unique solution is _uniq such that _uniq =:
		
    _uniq := (^* )^-1^* .

		 If there is no exact solution, the system is called an overdetermined inconsistent system. When a system is inconsistent and  is full rank, we often seek to minimize the sum of squared residuals, i.e. to find the ordinary least squares solution
		
    := (^* )^-1^* .

		The residual can be written as =  -. Note that ^* = 0, which can be easily seen by substituting =  + into (<ref>). For simplicity, we will refer to the matrix  of a linear system as consistent or inconsistent when the system itself is consistent or inconsistent.
		
		If the matrix  in the linear system = is rank deficient, then there are infinitely many solutions to the system regardless the size of m and n. In this case, we again want the least norm solution in the underdetermined case and the “least-norm least-squares" solution in the overdetermined case,
		
    :=
    ^*(^*)^† ,    if m < n
    
            (^* )^†^* ,    if m > n,

where (·)^† is the pseudo-inverse. General solutions to the linear system can be written as =  + where = 0—note that =  =  +  =. Similar to the full rank case, when the low-rank system is inconsistent, we can write =  +, again where ^* = 0. 

		

 §.§ Randomized Kaczmarz and its Extension

		
			The Kaczmarz Algorithm <cit.> solves a linear system = by cycling through rows of  and projects the estimate onto the solution space given by the chosen row. 
			It was initially proposed by Kaczmarz <cit.> and has recently regained interest in the setting of computer tomography where it is known as the Algebraic Reconstruction Technique <cit.>. 
			The randomized variant of the Kaczmarz method introduced by Strohmer and Vershynin <cit.> was proven to converge linearly in expectation for consistent systems.
			Formally, given  and  of (<ref>), RK chooses row i ∈{1, 2, ... m } of  with probability ^2_2/^2_F, 
			and projects the previous estimate onto that row with the update
			
    :=  + (_i - )/^2_2 ()^*.

			Needell <cit.> later studied the inconsistent case and showed that RK does not converge to the least squares solution for inconsistent systems, but rather converges linearly to some convergence radius of the solution.
			To remedy this, Zouzias and Freris <cit.> proposed the Randomized Extended Kaczmarz (REK) algorithm to solve linear systems in all settings. 
						For REK, row i∈{1, 2, ... m } and column j ∈{1, ...n } of  are chosen at random with probability 			
			
    (row = i) = ^2_2/^2_F ,   (column = j) = ^2_2/^2_F,

			and starting from _0=0 and _0=, every iteration computes 
			
    :=  + (^i -^i_t-  )/^2_2 ()^*,   _t := _t-1 - ⟨, _t-1⟩/_2^2.

			REK finds the optimal solution in all linear system settings. 
			In the consistent setting, it behaves as RK.
			In the overdetermined inconsistent setting,  estimates the residual vector  and allows  to converge to the true least squares solution of the system. 
			REK was shown to converge linearly in expectation to the least-squares solution by Zouzias and Freris <cit.>.
			



		

 §.§ Randomized Gauss-Seidel and its Extension

		The Gauss-Seidel method was originally published by Seidel but it was later discovered that Gauss had studied this method in a letter to his student <cit.>. 
		Instead of relying on rows of a matrix, the Gauss-Seidel method relies on columns of .
		The randomized variant was studied by Leventhal and Lewis <cit.> shortly after RK was published. 
		The randomized variant (RGS) requires a column j to be chosen randomly with probability ^2_2/^2_F,
	and updates at every iteration
		
    :=  + ^*(- )/^2_2_(j),

		where _(j) is the j^th basis vector (a vector with 1 in the j^th position and 0 elsewhere). Leventhal and Lewis <cit.> showed that RGS converges linearly in expectation when  is overdetermined.
		However, it fails to find the least norm solution for an underdetermined linear system <cit.>. 
				The Randomized Extended Gauss-Seidel (REGS) resolves this problem, much like REK did for RK in the case of overdetermined systems. 
		The method chooses a random row and column of  exactly as in REK, and then updates at every iteration
		 
    :=  + ^*(- )/^2_2_(j) ,
    _i    := _n - ()^*/_2^2 ,
    _t    := _i(_t-1 +  - ),

		and at any fixed time t, outputs - _t as the estimated solution to =. Here, _n denotes the n × n identity matrix. This extension works for all variations of linear systems and was proven to converge linearly in expectation by Ma et al. <cit.>.


				The RK and RGS methods along with their extensions are extensively studied and compared in <cit.>.
		Table <ref> summarizes the convergence properties of each of the randomized methods and their extensions. 
		
				
				



		In this paper, we focus on using combinations of RK and REK but also discuss RGS and REGS for comparison.
		We choose to focus on RK and REK because their updates consist only of scalar operations and inner products as opposed to REGS which requires an outer product. The methods proposed are easily  extendable to RGS and REGS. 

		
		




§ VARIATIONS OF FACTORED LINEAR SYSTEMS


Our proposed methods rely on interleaving solution estimates to subsystem (<ref>) and subsystem (<ref>).  
	Because the convergence of RK, RGS, REK, and REGS are heavily dependent on the number of rows and columns in the linear system, it is important to discuss how the settings of (<ref>) and  (<ref>) are determined by .
	In this section, we will discuss when we can expect our methods to solve the full system. 


	
	For simplicity in notation, we will denote  , , and  as the “optimal” solution of (<ref>), (<ref>), and (<ref>) respectively, as summarized in Table <ref>.
	By “optimal" solution for (<ref>) and (<ref>), we mean the unique, least norm, or the least squares solution, depending on the type of system (overdetermined consistent, underdetermined, overdetermined inconsistent). Since we assume that  may be low-rank,  is going to be the least norm solution as described in (<ref>).
	Table <ref> presents such a summary depending on the size of k with respect to m and n. 

	

 We spend the rest of this section justifying the shading in Table <ref>. 	For this, we split Table <ref> into three scenarios : (S1)  is overdetermined and consistent, (S2)  is underdetermined, (S3a and S3b)  is overdetermined and inconsistent. It should be noted that in Scenarios S1 and S2, the overdetermined-ness or underdetermined-ness of each subsystem follows immediately from sizes of m, n, and k and the assumption that the subsystems are full rank. We use the over/underdetermined-ness of each subsystem to show = for Scenario S1 and ≠ for Scenario S2. In Scenario S3, a little more work needs to be done to conclude the consistency of each subsystem. For Scenario S3a and S3b, we first investigate how inconsistency in (<ref>) affects the consistency of (<ref>) and (<ref>), then show that for Scenario S3a, ≠ and for Scenario S3b =.  This section provides the intuition on when we should expect our methods (or similar ones based on interleaving solutions to the subsystems) to work. However, one may also skip ahead to the next section where formally we present our algorithm and main results.
	
	
	
	
		
  * Scenario S1: U overdetermined, consistent.
			When  is overdetermined and consistent, we find that solving (<ref>) and (<ref>) gives us the optimal solution of (<ref>). 
			Indeed, in the case where  is overdetermined and consistent, we have:
			
    = (^*)^-1^* (^*)^-1^* 
       = (^*)^-1^* (^*)^-1^* 
        = (^*)^-1^* 
        =

			In the case where  is underdetermined, we have:
			
    = ^*(^*)^-1 (^*)^-1^* 
       = ^*(^*)^-1 (^*)^-1^* 
    
    						            =

				
			Since  is possibly low-rank, we still need to argue that this implies that =, i.e.  is indeed the least norm solution to the full system. Suppose towards a contradiction that  is the least norm solution of the full system but not subsystem (<ref>); in other words assume that =  + where = 0 and  is nontrivial. Multiplying both sides by , we see that
			since  has a nontrivial component  such that = 0, it  cannot be the least norm solution to the full system as assumed, reaching a contradiction. 
		Therefore, in the consistent case when  is overdetermined, we have that =, and may hope that our proposed methods will be able to solve the full system (<ref>) utilizing the subsystems (<ref>) and (<ref>).
		
  * Scenario S2: U underdetermined.
			When  is underdetermined, solving (<ref>) and (<ref>) for their optimal solutions does not guarantee the optimal solution of the full system. 
			Intuitively, (<ref>) has infinitely many solutions and = _LN≠.
			Mathematically, investigating , we find that 
			
    = ^*(^*)^-1^*(^*)^-1
       = ^*(^*)^-1^*(^*)^-1
       = ^*(^*)^-1
       ≠  if   underdetermined 
        = (^*)^-1^* (^*(^*)^-1 - _V)
       = (^*)^-1^* ^*(^*)^-1
       ≠  if   overdetermined, inconsistent,

			where we rewrite =  + _V for _V ∈ null(^*) since subsystem (<ref>) may be inconsistent. Note that if subsystem (<ref>) is consistent then we simply have _V = 0 and the above calculation still carries through.
Therefore, we do not expect our proposed methods to succeed when  is underdetermined.  Fortunately, this case seems to be of little practical interest, since factoring an underdetermined system does not typically save any computation. 
		
  * Scenario S3: X inconsistent. 
			Before we discuss whether it's possible to recover the optimal solution to the full system, we must first discuss what  being inconsistent implies about the subsystems (<ref>) and (<ref>).
			In particular, one needs to determine whether inconsistency in the full system creates inconsistencies in the individual subsystems. 
			If  is inconsistent then we have +  = where  is the optimal solution of (<ref>) and ^*  = 0.
			Now, consider decomposing = _1 + _2 where ^* _1 = 0, ^* _2 ≠0, and ^*^*_2 = 0. 
			Notice that ^*  = ^*^*(_1 + _2) = ^*^*_1 + ^*^*_2 = 0, as desired. 
			We want to decompose the full system +  = into two subsystems.
			Following a similar thought process as before, we choose to decompose our full system into the following: 
			
    + _1 + _2= 
        = .


			Clearly, (<ref>) is inconsistent since ^*_1=0 and ^* _2 ≠0. 
			
			Because  must be overdetermined for  to be inconsistent, = (^* )^-1^*( - _1 - _2) = (^* )^-1^*( -  _2) is the least squares solution to (<ref>).
			
				
				
  * Case S3a:  overdetermined. 
				Note that the second subsystem (<ref>) is possibly inconsistent (since there may be a component  of in the null space of ^*).
					Writing =  + _V such that _V ∈ null(^*), we have
				    
    = (^*)^-1^* ((^* )^-1^* ( - _1 - _2) - _V )
       = (^*)^-1^* (^* )^-1^* ( - _2) 
       =  - (^*)^-1^* (^* )^-1^*_2 
       ≠.

				   Similar to Scenario S2, if subsystem (<ref>) is indeed consistent then _V = 0 and the above calculation still carries through.
				   Therefore, in this case, we do not expect to find the optimal solution to (<ref>). 
				   
				   
  * Case S3b:  underdetermined. 
					In this case, _2 = 0 and solving (<ref>) and (<ref>) obtains the optimal solution to the full system since
					
    = ^*(^*)^-1 (^* )^-1^*( - _1) 
        = ^*(^*)^-1 (^* )^-1^* 
        = ^*(^*)^-1 (^* )^-1^* 
        =^*(^*)^-1
        =
				
					Following the same argument as in Scenario S1 when  is underdetermined, we reach the conclusion that =. Thus, in this case our methods have the potential to solve the full system.

				
					
	

	These three scenarios fully explain the shading in Table <ref>. 
	The focus of the remainder of this paper will be the case in which k < m,n (i.e. left column of Table <ref>) since, as mentioned, it is  practically the most relevant setting.
		
			




§ METHODS AND MAIN RESULTS
 
	
	
	Our approach intertwines two iterative methods to solve subsystem (<ref>) followed by subsystem (<ref>). 
	For the consistent setting, we propose Algorithm <ref> which uses an iterate of RK on (<ref>) intertwined with an iterate of RK to solve (<ref>). 
	For the inconsistent setting, we propose using REK to solve subsystem (<ref>) followed by RK to solve subsystem (<ref>) as shown in Algorithm <ref>. We view the latter method as a more practical approach, and the former as interesting from a theoretical point of view.
	Recall that ^p is the p^th element in the vector .
	Standard stopping criteria include terminating when the difference in the iterates is small or when the residual is less than a predetermined tolerance. To avoid adding complexity to the algorithm, the residual should be computed and checked approximately every m iterations.
		We propose an approach that interlaces solving subsystems (<ref>) and (<ref>); this has a couple of advantages over solving each subsystem separately.
	First, if we are given some tolerance ϵ that we allow on the full system, it is unclear when we should stop the iterates of the first subsystem to obtain such an error — if solving the first subsystem is terminated prematurely, the error may propagate through iterates when solving the second subsystem.
	Second, the interlacing allows for opportunities to implement these algorithms in parallel. 
	We leave the specifics of such an implementation as future work as it is outside the scope of this paper.

	
	
	
	
	
	
	
	

 §.§ Main result

	
	Our main result shows that Algorithm <ref> and Algorithm <ref> converge linearly to the desired solution. The convergence rate, as expected, is a function of the conditioning of the subsystems, and hence we introduce the following notation. Here and throughout, for any matrix  we write 
	
    α_A     :=  1 - σ_min^2()/^2_F , 	
        := σ_max^2()/σ_min^2() ,
 
	   where σ_min^2() is the smallest non-zero singular value of , and  is the squared condition number of .
	  Recall that the optimal solution to a system is either the least-norm, unique, or least-squares solution depending on whether the system is underdetermined, overdetermined consistent, or overdetermined inconsistent, respectively.
	
		
		
		 Let  be low rank, = such that ∈ℂ^m × k and ∈ℂ^k × n are full rank, and the systems =, =, and = have optimal solutions , , and  respectively, and , , κ^2_U are as defined in (<ref>) and (<ref>). Setting _0 = 0 and assuming k < m,n, we have
			
				
  * if  = is consistent, then = and Algorithm <ref> converges with expected error
					
    - ^2 ≤^t ^2 +  (1-γ_1)^-1α_max^t^2/^2_F,     if ≠
    ^t ^2 + t α_max^t ^2/^2_F ,       else

					where α_max = max{, } and γ_1 = min{/, /}.
				
  * if = is inconsistent, then  = and Algorithm <ref> converges with expected error
					
    - ^2 ≤^t ^2 +(1-γ_2)^-1α̃_max^t-1(1 + 2κ_U^2)^2/^2_F,      if √()≠
    ^t ^2 + t α̃_max^t-1(1 + 2κ_U^2)^2/^2_F,                 else

					where α̃_max = max{√(), } and γ_2 = min{√()/, /√()}.
			
			
		
		
		Remarks. 

		1. Theorem <ref>(a) also applies to the setting in which  is overdetermined, consistent, and n < k < m. In the proof of Lemma <ref>, one must simply note that the bound (<ref>) still holds for this setting and all other steps in the proof analogously follow.

	2. Empirically, our experiments in the next section suggest that  and  can be substantially better conditioned than . 

 	3. Algorithm <ref> is interesting to discuss from a theoretical standpoint but in applications Algorithm <ref> is more practical as linear systems are typically inconsistent. Algorithm <ref> can be utilized in applications if error in the solution is tolerable. In particular, if = is inconsistent then Algorithm <ref> will converge in expectation to some convergence horizon. This can be see by replacing the use of Proposition <ref> in the bound (<ref>) with the convergence bound of RK on inconsistent linear systems found in Theorem 2.1 of <cit.>.

 	4. While not the main focus of this paper, we briefly note here that for matrices large enough that they cannot be stored entirely in memory, there is an additional cost that must be paid in terms of moving data between the disk and RAM. In a truly large-scale implementation, the RK-RK algorithm might be more scalable than the REK-RK algorithm since RK only accesses random rows of both  and , which is efficient if both matrices are stored in row major form, but REK accesses both random rows and columns, and hence storing in either row major or column major format will be slow for one of the two operations.



 §.§ Supporting results


		To prepare for the proof of the above theorem (the central theoretical result of the paper), we state a few supporting results which will help simplify the presentation of the proof. We begin by stating known results on the convergence of RK and REK on linear systems. 
		Let _b denote the expected value taken over the choice of rows in  and _x the expected value taken over the choice of rows in  and when necessary the choice of columns in .
		Also, let  denote the full expected value (over all random variables and iterations) and ^t-1 be the expectation conditional on the first t-1 iterations.
		
		(<cit.>) Given a consistent linear system =, the Randomized Kaczmarz algorithm, with initialization _0 = 0, as described in Section <ref> converges to the optimal solution  with expected error
			
    _t - ^2 ≤^t^2,

			where  is as defined in (<ref>).
			
		
		
		(<cit.>) Given a linear system =, the Randomized Extended Kaczmarz algorithm, with initialization _0 = 0, as described in Section <ref> converges to the optimal solution  with expected error
			
    _t - ^2 ≤^⌊ t/2 ⌋ (1 + 2 ) ^2,

			where  is as defined in (<ref>) and  is as defined in (<ref>).
			
		
		
		The proof of Theorem <ref> builds directly on two useful lemmas. Lemma <ref> addresses the impact of intertwining the algorithms. In particular, it shows useful relationships involving , the RK update solving the linear system = at the t^th iteration (with  as the previous estimate), and our update . Lemma <ref> states that conditional on the first t-1 iterations, we can split the norm squared error - ^2 into two terms relating to the error from solving subsystem (<ref>) and the error from solving subsystem (<ref>). To complete the proof of Theorem <ref>, we bound the error from solving (<ref>) depending on whether we use RK (as in Algorithm <ref>) or REK (as in Algorithm <ref>) then apply the law of iterated expectations to bound the error from solving (<ref>). We now state the aforementioned lemmas, and then formally prove the theorem.
		
		 Let =  + ^p - /^2  ()^*. In Algorithm <ref> and Algorithm <ref> we have that:
		
		
  * _b^t-1⟨ - ,  - ⟩ = 0, 
		
  * - ^2 =  - ^2 -  - ^2.
		
			
		 
		In words, part (a) states that the difference between an RK iterate solving the exact linear system = and our RK iterate (which solves the linear system resulting from intertwining =), is orthogonal to -. This will come in handy in Lemma <ref>. Part (b) is a Pythagoras-style statement, which follows from well-known orthogonality properties of RK updates, included here for simplicity and completeness.
		

		To prove statement (b), we note that ( - ) is parallel to  and ( - ) is perpendicular to  since ( - ) = ( + ^p - /^2  ()^* - ) =  + ^p -  - ^p = 0. We apply the Pythagorean Theorem to obtain the desired result.

		We prove statement (a) by direct substitution and expansion, as follows:
		
    _b^t-1   ⟨ - ,  - ⟩ 
    			= _b^t-1⟨^p - ^p/^2  ()^*,  -  + ^p- /^2 ()^* ⟩
       (i)=_b^t-1⟨^p - ^p/^2  ()^*, ^p- /^2 ()^* ⟩  +  _b^t-1⟨^p - ^p/^2 ()^*,  - ⟩
       (ii)=_b^t-1 (^p - ^p)( ^p-  )  /^2   +  ⟨_b^t-1^p - ^p/^2  ()^*,  - ⟩
       (iii)=∑_p  (^p - ^p)( ^p-  )  /^2 ^2/^2_F  +  ⟨∑_p   (^p - ^p) ()^* /^2 ^2/^2_F,  - ⟩
       =  ( - )^*(-  )  /^2_F  +  ⟨^*( - ) /^2_F,  - ⟩
       (iv)=⟨ - /^2_F  , (- ) ⟩ + ⟨ - /^2_F   , (  - )⟩
       = 0.

			Step (i) follows from linearity of inner products, step (ii) simplifies the inner product of two parallel vectors, and step (iii) computes the expectation over all possible choices of rows of . In step (iv), we use the fact that for k < m,n, subsystem (<ref>) is always consistent (since  is underdetermined) to make the substitution =. 

				
 		
		 In Algorithm <ref> and Algorithm <ref>, we can bound the expected norm squared error of - as
			
    ^t-1 - ^2  ≤ - ^2 + _x^t-1 - ^2/^2_F .
 
			
		
We investigate the expectation of the norm squared error of - conditional on the first t-1 iterations and over the choice of rows of . We keep _x^t-1 in our bound as this expectation will depend on whether Algorithm <ref> or Algorithm <ref> is being used. 
		
		
		
    ^t-1 - ^2 
    				   =^t-1 -  +  - ^2 
       =^t-1 - ^2  + ^t-1 - ^2 + 2 ^t-1⟨ - ,  - ⟩
       (iii)=^t-1 - ^2  + ^t-1 - ^2 
       (iv)=^t-1 - ^2 - ^t-1 - ^2 + ^t-1 - ^2 
       (v)= - ^2 -  ^t-1 (^p - )  /^2  ()^* ^2  + ^t-1(^p - ^p)/^2   ()^* ^2  
        =  - ^2 -  ^t-1[  |- |^2  /^2 ]  + ^t-1[ |^p - ^p|^2/^2 ].

				Steps (iii) and (iv) are applications of Lemma <ref>(a) and Lemma <ref>(b) respectively, and step (v) follows from the definition of each term and simplification using the fact that =.
	 
	 Now, we evaluate the conditional expectation on the choices of rows of  to complete the proof:
	 
    ^t-1 - ^2 
    			   (vi)= - ^2 -  _b^t-1[  |- |^2  /^2 ]  + _x^t-1_b^t-1[ |^p - ^p |^2/^2 ] 
       = - ^2 -  ∑_p = 1^k  |- |^2  /^2 ^2/^2_F  + _x^t-1∑_p = 1^k |^p - ^p |^2/^2 ^2/^2_F
       = - ^2 -  - ^2  /^2_F  + _x^t-1[  - ^2/^2_F ] 
       (vii)≤ - ^2 - σ_min^2()  -  ^2  /^2_F  + _x^t-1[ - ^2/^2_F ] 
        =  - ^2 + _x^t-1[  - ^2/^2_F ] .

		In step (vi), we use iterated expectations to split the expected value ^t-1 = _x^t-1_b^t-1. Step (vii) uses the fact that ( - ) ^2 ≥σ_min^2()  - ^2 since - are in the row span of  for all t. We simplify and obtain the desired bound.
	 	 
	


 §.§ Proof of main result


	We now have all the ingredients we need to prove Theorem <ref>, which we now proceed to below.
	
	
The fact that = was already argued in scenarios S1 and S3(b) in the previous section, so we do not reproduce its argument here. Given this fact, to prove Theorem <ref>, we only need to invoke the statement of Lemma <ref> and bound the term _x^t-1 - ^2 using Proposition <ref> or Proposition <ref> depending on whether we are using Algorithm <ref> or Algorithm <ref>, respectively. 
		 
			
  * For Algorithm <ref>, plugging Proposition <ref> into the statement of Lemma <ref> yields
				
    ^t-1 - ^2 
    					
    					   ≤ - ^2 + ^t ^2/^2_F .

					Let α_max = max{α_U, α_V }, γ_1 = min{α_U/α_V, α_V/α_U}, and note that γ_1 < 1 if α_U ≠α_V. Taking expectations over the randomness from the first t-1 iterations and using the Law of Iterated Expectation, we have
					
    - ^2    ≤^t ^2 + ^2/^2_F ∑_h=0^t-1^t-h^h 
       =^t ^2 + ^2/^2_F ∑_h=0^t-1^t-1-h^h 
       (i)=^t ^2 + α_max^t-1^2/^2_F ∑_h=0^t-1γ_1^h
       (ii)≤^t ^2 + α_max^t^2/^2_F 1/1-γ_1.

				where step (i) uses the fact that the summation is symmetric with respect to α_U and α_V. Step (ii) uses the fact that α_max^t-1 = α_max^t if α_max = and α_max^t-1 < α_max^t if α_max =. Now, if α_U = α_V = α_max then we have
					
    - ^2    ≤^t ^2 + ^2/^2_F ∑_h=0^t-1α_max^t-hα_max^h 
       (ii)=^t ^2 + tα_max^t ^2/^2_F ,

					where the second term of (ii) approaches 0 as t →∞ since α_max< 1.	
			
  * For Algorithm <ref>, plugging Proposition <ref> into the statement of Lemma <ref> yields
				
    ^t-1 - ^2 
    					
    					   ≤ - ^2 + (1+2)^⌊ t/2 ⌋^2/^2_F .

					Taking expectations over the remaining randomness, and using the fact that ^⌊t-h/2⌋≤^t-h-1/2 since < 1, we have
					
    - ^2  
    					   ≤^t ^2 + (1+2) ^2/^2_F ∑_h=0^t-1^⌊t-h/2⌋^h 
       ≤^t ^2 + (1+2) ^2/^2_F ∑_h=0^t-1^t-h-1/2^h 
        =  ^t ^2 + (1+2) ^2/^2_F ∑_h=0^t-1√()^t-1-h^h .

					Using the same techniques as in (a), let α̃_max = max{, √()} and γ_2 = min{/√() , √()/} and again noting γ_2 < 1 when ≠√(), we can write:
					
    - ^2  
    					   ≤^t ^2 + (1+2) ^2/^2_F ∑_h=0^t-1√()^t-1-h^h 
        =  ^t ^2 + (1+2) α̃_max^t-1^2/^2_F ∑_h=0^t-1γ_2^h 
       ≤^t ^2 + α̃_max^t-1 (1+2) ^2/^2_F 1/1-γ_2  .

				When = √() = α̃_max, 
				
    - ^2  
    					   ≤^t ^2 + (1+2) ^2/^2_F ∑_h=0^t-1α̃_max^t-1-hα̃_max^h 
       (i)=^t ^2 + tα̃_max^t-1(1+2)  ^2/^2_F ,

				where the second term of (i) goes to 0 as t goes to infinity.
				
		This concludes the proof of the theorem.
	





§ EXPERIMENTS

In this section we discuss experiments done on both simulated and real data using different algorithms in different settings. 
The naming convention for the remainder of the paper will be to refer to ALG1-ALG2 as an interlaced algorithm where ALG1 is the algorithm iterate used to solve subsystem (<ref>) and ALG2 is the algorithm used to solve subsystem (<ref>). When an algorithm's name is used alone, we imply applying the algorithm on the full system (<ref>).





In Figure <ref> we show our first set of experiments. Entries of , , and  are drawn from a standard Gaussian distribution. We set = and = if  is consistent and =  + where ∈ null (^*) (computed in Matlab using  function) if  is inconsistent. In this first set of experiments, m,n,k ∈{ 100,   150,   200 } depending on the desired size of k with respect to the over or underdetermined-ness of . For example, if k < m,n and  is overdetermined then k = 100, m = 200, and n = 150. The plots show iteration vs ℓ_2-error, - ^2, of each method averaged over 40 runs and allowing each algorithm to run 7 × 10^4 iterations. The layout of Figure <ref> is exactly as in Table <ref>. For each row, we have a different setting for  and for each column, we vary the size of k depending on the size of . Looking at the overall trends, we see that when k < m,n and when  is overdetermined, consistent and n < k < m, there is a method that obtains the optimal solution for the system. These results align with the expectations set in Table <ref>. Looking at each individual subplot, we also find what one would expect according to Table <ref>. In other words, if  or  is in one of the settings where RK or RGS are expected to fail then RK-RK or RGS-RGS fail as well.

When  is overdetermined, inconsistent and k < m,n we have that  is underdetermined. In this case, we don't need to interlace iterates of REK and REK together. To work on an underdetermined system, using RK is enough to find the optimal solution of that subsystem. This motivated interlacing iterates of RK with REK. Figure <ref> has the same set up as discussed in the previous experiment with the exception of using larger random matrices with : 1200 × 750 and k = 500. In Figure <ref> we plot iteration vs ℓ_2-error and in Figure <ref> we plot FLOPS vs ℓ_2-error. The errors are averaged over 40 runs, with shaded regions representing error within two standard deviations from the mean. Note that Algorithm <ref> performs excellently in practice, better than our theoretical upper bound as computed in Theorem <ref>. In this experiment, we see that REK-RK and REK-REK perform comparably in error and that REK-RK is more efficient in FLOPS.



In addition to simulated experiments, we also show the usefulness of these algorithms on real world data sets on wine quality, bike rental data, and Yelp reviews. In all following experiments, we plot the average ℓ_2-error at the t^th iteration over 40 runs and shaded regions representing the ℓ_2-error within two standard deviations. In addition to empirical performance, we also plot the theoretical convergence bound derived in Theorem <ref> (labeled “BND" in the legends). From these experiments, it is clear that the algorithms perform even better in practice than the worst-case theoretical upper bound. The data sets on wine quality and bike rental data are obtained from the UCI Machine Learning Repository <cit.>. The wine data set is a sample of m=1599 red wines with n = 11 physio-chemical properties of each wine. We choose k = 5 and compute  and  using Matlab's  function for nonnegative matrix factorization (recall the motivations from the first section). Figure <ref> shows the results from this experiment. The conditioning of , , and  are = 2.46 × 10^3, = 25.96, and = 4.20 respectively. We plot the ℓ_2-error averaged over 40 runs. Since  has such a large condition number, this impacts the convergence of REK on  negatively as shown by the seemingly horizontal line (the error is actually decreasing, but incredibly slowly). We also see that REK-RK and REK-REK are working comparably and significantly faster than REK alone. This can be explained by the better conditioning on  and .



The bike data set contains hourly counts of rental bikes in a bike share system. The data sets contains date as well as weather and seasonal data. There are m = 17379 samples and n = 9 attributes per sample. We choose k = 8 and compute  and  in the same way as with the wine data set. Figure <ref> shows the results from this experiment. The conditioning of , , and  are = 94.27, = 54.91, and = 2.99 respectively. Similar to Figure <ref>, we see that the convergence of REK suffers from the poorly conditioned matrix . We also see again that REK-REK and REK-RK behave similarly and outperform REK.

To show the advantage of our algorithms on large systems, we create extremely large standard Gaussian matrices : 10^6 × 10^3 and : 10^3 × 10^4. These matrices are so large that the matrix product  cannot be computed in Matlab due to memory constraints. These results are shown in Figure <ref>. We see that without needing to do the matrix computation, we are still able to find the solution to the linear system =.



Lastly, we present the performance of our methods on a large real world data set. We use the Yelp challenge data set <cit.>. In our setting, we let : 10^5 × 10^4 be a document term frequency matrix where each row represents a Yelp review and each column represents a word feature. The elements of  contain the frequency at which the word is used in the review. We only use a subset of the amount of data available due Matlab memory constraints. Here,  is a vector that represents the number of stars a review received. We choose k = 5000. Figure <ref> shows the results from this experiment using REK, REK-REK, and REK-RK. The conditioning of , , and  are  = 127.3592, = 24.274, and = 19.096 respectively. In this large real world data set, we can again see the usefulness of our proposed methods when we are given =. 





These experiments complement and verify our theoretical findings. In settings which we expect to fail to obtain the least squares or least norm solutions, our experiments show that they do indeed fail. Additionally, where we expect that the optimal solution is obtainable, the experiments 
show the proposed methods can obtain such solutions and in many instances outperform the original algorithm on the full system. We see that empirically, subsystems are better conditioned than full systems, thus explaining their better performance.






§ CONCLUSION




We have proposed two methods interlacing Kaczmarz updates to solve factored systems.  For large-scale applications in which the system is stored in factored form for efficiency or the factorization arises naturally, our methods allow one to solve the system without the need to perform the large-scale matrix product first.  Our main result proves that our methods provide linear convergence in expectation to the (least-squares or least-norm) solution of (overdetermined or underdetermined) linear systems.  Our experiments support these results, and show that our methods provide significant computational advantages for factored systems.  The interlaced structure of our methods suggests they can be implemented in parallel which would lead to even further computational gains. We leave such details for future work.  Additional future work includes the design and analysis of methods that converge to the solution in the settings not covered in this paper, i.e. the gray cells of Table <ref>.  Although its practical implications are not immediately clear to us, these may still be of theoretical interest.







§ ACKNOWLEDGMENTS

Needell was partially supported by NSF CAREER grant #1348721, and the Alfred P. Sloan Fellowship. Ma was supported in part by NSF CAREER grant #1348721, the CSRC Intellisis Fellowship, and the Edison International Scholarship. The authors would like to thank Wutao Si for pointing out a flaw in the original proof of Theorem 1.  In addition, the authors also thank the Institute of Pure and Applied Mathematics (IPAM) where this collaboration started. 
abbrvnat

