Iterative methods for solving factorized linear systems

TL;DR

This paper introduces a new iterative algorithm for solving large linear systems stored in factorized form, achieving faster convergence by leveraging the matrix factorization without explicitly computing the full matrix.

Contribution

It proposes a novel randomized Kaczmarz variant that exploits matrix factorization, providing theoretical convergence guarantees and practical acceleration.

Findings

Exponential convergence rate proven for the proposed method.

Significant acceleration observed in experiments.

Method effectively handles large low-rank datasets.

Abstract

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 X that is stored in a factorized form, X = UV; 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 X. We prove an exponential convergence rate and supplement our theoretical guarantees with experimental evidence demonstrating that the factored variant yields significant acceleration in convergence.

Tables

Table 1. Table 1 . Summary of convergence properties of randomized methods under all settings.

Method	Overdetermined, consistent : convergence to ${𝜷}_{\text{uniq}}$?	Overdetermined, inconsistent : convergence to ${𝜷}_{𝑳𝑺}$?	Underdetermined : convergence to ${𝜷}_{𝑳𝑵}$?
RK	Yes [25]	No [22]	Yes [19]
REK	Yes [30]	Yes [30]	Yes [19]
RGS	Yes [16]	Yes [16]	No [19]
REGS	Yes [19]	Yes [19]	Yes [19]

Table 2. Table 2 . Summary of notation for linear systems discussed and their solutions

Linear System	Optimal Solution
$𝑿𝜷=𝒚$ (1)	${𝜷}_{\star }$
$𝑼𝒙=𝒚$ (2)	${𝒙}_{\star }$
$𝑽𝒃=𝒙$ (3)	${𝒃}_{\star }$

Table 3. Table 3 . Summary of types of matrices 𝑼 𝑼 {\boldsymbol{U}} and 𝑽 𝑽 {\boldsymbol{V}} for given m , n , 𝑚 𝑛 m,n, and k 𝑘 k relations. The first column indicates the setting where k 𝑘 k is less than both m 𝑚 m and n 𝑛 n , the second column when k 𝑘 k is between m 𝑚 m and n 𝑛 n , and the third when k 𝑘 k is greater than both m 𝑚 m and n 𝑛 n . The cells in white indicate when our proposed methods will converge to the solution of the full system, gray cells indicate when our methods will not. Recall that 𝑿 ∈ ℂ m × n 𝑿 superscript ℂ 𝑚 𝑛 {\boldsymbol{X}}\in\mathbb{C}^{m\times n} , 𝑼 ∈ ℂ m × k 𝑼 superscript ℂ 𝑚 𝑘 {\boldsymbol{U}}\in\mathbb{C}^{m\times k} and 𝑽 ∈ ℂ k × n 𝑽 superscript ℂ 𝑘 𝑛 {\boldsymbol{V}}\in\mathbb{C}^{k\times n} . We denote with “(In)con” systems that can be either consistent or inconsistent. Arguably, the k < min ⁡ { m , n } 𝑘 𝑚 𝑛 k<\min\{m,n\} setting is most practically relevant, and in this case our methods do indeed recover the optimal solution.

$𝑿$					$k>\max \{m,n\}$
Underdetermined	$𝑼$ = Over, Consis. $𝑽$ = Under	(S1)	$𝑼$ = Under $𝑽$ = Under	(S2)	$𝑼$ = Under $𝑽$ = Over, (In)con.	(S2)
Overdetermined Consis.	$𝑼$ = Over, Consis. $𝑽$ = Under	(S1)	$𝑼$ = Over, Consis. $𝑽$ = Over, Consis.	(S1)	$𝑼$ = Under $𝑽$ = Over, (In)con.	(S2)
Overdetermined Inonsis.	$𝑼$ = Over, Incon. $𝑽$ = Under	(S3b)	$𝑼$ = Over, Incon. $𝑽$ = Over, (In)con.	(S3a)	$𝑼$ = Under $𝑽$ = Over, (In)con.	(S2)

Equations

$$\mathbf{U}\mathbf{V}\mathbf{\beta }=\mathbf{y},$$

$$\mathbf{U}\mathbf{x}$$

$$\mathbf{V}\mathbf{\beta }$$

$${\mathbf{\beta }}_{\mathbf{LN}}:={\mathbf{X}}^{\ast }(\mathbf{X}{\mathbf{X}}^{\ast }{)}^{-1}\mathbf{y}.$$

$${\mathbf{\beta }}_{\text{uniq}}:=({\mathbf{X}}^{\ast }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\ast }\mathbf{y}.$$

$${\mathbf{\beta }}_{\mathbf{LS}}:=({\mathbf{X}}^{\ast }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\ast }\mathbf{y}.$$

$${\mathbf{\beta }}_{\mathbf{LN}}:=\{\begin{array}{ll}{\mathbf{X}}^{\ast }(\mathbf{X}{\mathbf{X}}^{\ast }{)}^{†}\mathbf{y},& \text{if m<n }\\ ({\mathbf{X}}^{\ast }\mathbf{X}{)}^{†}{\mathbf{X}}^{\ast }\mathbf{y},& \text{if m>n },\end{array}$$

$${\mathbf{\beta }}_{\mathbf{t}}:={\mathbf{\beta }}_{\mathbf{t}-1}+\frac{({\mathbf{y}}_i-{\mathbf{X}}^i{\mathbf{\beta }}_{\mathbf{t}-1})}{{\Vert {\mathbf{X}}^i\Vert }_2^2}({\mathbf{X}}^i{)}^{\ast }.$$

$$\mathbf{P}(row=i)=\frac{{\Vert {\mathbf{X}}^i\Vert }_2^2}{{\Vert \mathbf{X}\Vert }_F^2},\mathbf{P}(column=j)=\frac{{\Vert {\mathbf{X}}_{(j)}\Vert }_2^2}{{\Vert \mathbf{X}\Vert }_F^2},$$

$${\mathbf{\beta }}_{\mathbf{t}}:={\mathbf{\beta }}_{\mathbf{t}-1}+\frac{({\mathbf{y}}^i-{\mathbf{z}}_t^i-{\mathbf{X}}^i{\mathbf{\beta }}_{\mathbf{t}-1})}{\Vert {\mathbf{X}}^i{\Vert }_2^2}({\mathbf{X}}^i{)}^{\ast },{\mathbf{z}}_t:={\mathbf{z}}_{t-1}-\frac{⟨{\mathbf{X}}_{(j)},{\mathbf{z}}_{t-1}⟩}{\Vert {\mathbf{X}}_{(j)}{\Vert }_2^2}{\mathbf{X}}_{(j)}.$$

$${\mathbf{\beta }}_{\mathbf{t}}:={\mathbf{\beta }}_{\mathbf{t}-1}+\frac{{{\mathbf{X}}_{(j)}}^{\ast }(\mathbf{y}-\mathbf{X}{\mathbf{\beta }}_{\mathbf{t}-1})}{\Vert {\mathbf{X}}_{(j)}{\Vert }_2^2}{\mathbf{e}}_{(j)},$$

$${\mathbf{\beta }}_{\mathbf{t}}$$

$${\mathbf{P}}_i$$

$${\mathbf{z}}_t$$

$${\mathbf{b}}_{\star }$$

$$=({\mathbf{V}}^{\ast }\mathbf{V}{)}^{-1}{\mathbf{V}}^{\ast }({\mathbf{U}}^{\ast }\mathbf{U}{)}^{-1}{\mathbf{U}}^{\ast }\mathbf{U}\mathbf{V}{\mathbf{\beta }}_{\star }$$

$$=({\mathbf{V}}^{\ast }\mathbf{V}{)}^{-1}{\mathbf{V}}^{\ast }\mathbf{V}{\mathbf{\beta }}_{\star }$$

$$={\mathbf{\beta }}_{\star }$$

$$\mathbf{V}{\mathbf{b}}_{\star }$$

$$=\mathbf{V}{\mathbf{V}}^{\ast }(\mathbf{V}{\mathbf{V}}^{\ast }{)}^{-1}({\mathbf{U}}^{\ast }\mathbf{U}{)}^{-1}{\mathbf{U}}^{\ast }\mathbf{U}\mathbf{V}{\mathbf{\beta }}_{\star }$$

$$=\mathbf{V}{\mathbf{\beta }}_{\star }$$

$$\mathbf{V}{\mathbf{b}}_{\star }$$

$$=\mathbf{V}{\mathbf{V}}^{\ast }(\mathbf{V}{\mathbf{V}}^{\ast }{)}^{-1}{\mathbf{U}}^{\ast }(\mathbf{U}{\mathbf{U}}^{\ast }{)}^{-1}\mathbf{U}\mathbf{V}{\mathbf{\beta }}_{\star }$$

$$={\mathbf{U}}^{\ast }(\mathbf{U}{\mathbf{U}}^{\ast }{)}^{-1}\mathbf{U}\mathbf{V}{\mathbf{\beta }}_{\star }$$

$$\ne \mathbf{V}{\mathbf{\beta }}_{\star }\text{ if V underdetermined }$$

$${\mathbf{b}}_{\star }$$

$$=({\mathbf{V}}^{\ast }\mathbf{V}{)}^{-1}{\mathbf{V}}^{\ast }{\mathbf{U}}^{\ast }(\mathbf{U}{\mathbf{U}}^{\ast }{)}^{-1}\mathbf{U}\mathbf{V}{\mathbf{\beta }}_{\star }$$

$$\ne {\mathbf{\beta }}_{\star }\text{ if V overdetermined, inconsistent, }$$

$$\mathbf{U}\mathbf{x}+{\mathbf{r}}_1+{\mathbf{r}}_2=\mathbf{y}$$

$$\mathbf{V}\mathbf{b}=\mathbf{x}.$$

$${\mathbf{b}}_{\star }$$

$$=({\mathbf{V}}^{\ast }\mathbf{V}{)}^{-1}{\mathbf{V}}^{\ast }({\mathbf{U}}^{\ast }\mathbf{U}{)}^{-1}{\mathbf{U}}^{\ast }(\mathbf{y}-{\mathbf{r}}_2)$$

$$={\mathbf{\beta }}_{\star }-({\mathbf{V}}^{\ast }\mathbf{V}{)}^{-1}{\mathbf{V}}^{\ast }({\mathbf{U}}^{\ast }\mathbf{U}{)}^{-1}{\mathbf{U}}^{\ast }{\mathbf{r}}_2$$

$$\ne {\mathbf{\beta }}_{\star }.$$

$$\mathbf{V}{\mathbf{b}}_{\star }$$

$$=\mathbf{V}{\mathbf{V}}^{\ast }(\mathbf{V}{\mathbf{V}}^{\ast }{)}^{-1}({\mathbf{U}}^{\ast }\mathbf{U}{)}^{-1}{\mathbf{U}}^{\ast }\mathbf{y}$$

$$=\mathbf{V}{\mathbf{V}}^{\ast }(\mathbf{V}{\mathbf{V}}^{\ast }{)}^{-1}({\mathbf{U}}^{\ast }\mathbf{U}{)}^{-1}{\mathbf{U}}^{\ast }\mathbf{U}\mathbf{V}{\mathbf{\beta }}_{\star }$$

$$=\mathbf{V}{\mathbf{V}}^{\ast }(\mathbf{V}{\mathbf{V}}^{\ast }{)}^{-1}\mathbf{V}{\mathbf{\beta }}_{\star }$$

$$=\mathbf{V}{\mathbf{\beta }}_{\star }$$

$${\alpha }_A$$

$${\kappa }_A^2$$

$$\mathbb{E}\Vert {\mathbf{b}}_t-{\mathbf{\beta }}_{\star }{\Vert }^2\le \{\begin{array}{l}{\alpha }_V^t\Vert {\mathbf{b}}_{\star }{\Vert }^2+(1-{\gamma }_1{)}^{-1}{\alpha }_{\text{max}}^t\frac{\Vert {\mathbf{x}}_{\star }{\Vert }^2}{\Vert \mathbf{V}{\Vert }_F^2},\text{ if }{\alpha }_U\ne {\alpha }_V\\ {\alpha }_V^t\Vert {\mathbf{b}}_{\star }{\Vert }^2+t{\alpha }_{\text{max}}^t\frac{\Vert {\mathbf{x}}_{\star }{\Vert }^2}{\Vert \mathbf{V}{\Vert }_F^2},\text{ else}\end{array}$$

Full text

1. Introduction

Recently, revived interest in stochastic iterative methods like the Kaczmarz [11, 25, 22, 23] and Gauss-Seidel [6, 19] 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 ${\boldsymbol{X}}{\boldsymbol{\beta}}={\boldsymbol{y}}$ with ${\boldsymbol{X}}\in\mathbb{C}^{m\times n}$. However, instead of the full system ${\boldsymbol{X}}$, we only have access to ${\boldsymbol{U}},{\boldsymbol{V}}$ such that ${\boldsymbol{X}}={\boldsymbol{U}}{\boldsymbol{V}}$. In this case, we want to solve the linear system:

[TABLE]

where ${\boldsymbol{U}}\in\mathbb{C}^{m\times k}$ and ${\boldsymbol{V}}\in\mathbb{C}^{k\times n}$. Instead of taking the product of ${\boldsymbol{U}}$ and ${\boldsymbol{V}}$, to form ${\boldsymbol{X}}$, which may not be desirable, we approach this problem using stochastic iterative methods to solve the individual subsystems

[TABLE]

in an alternating fashion. Note that ${\boldsymbol{\beta}}$ in (3) is the vector of unknowns that we want to solve for in (1) and ${\boldsymbol{y}}$ in (2) is the known right hand side vector of (1). If we substitute (3) into (2), we acquire the full linear system (1). We will often refer to (1) as the “full system” and (2) and (3) 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 ${\boldsymbol{x}}$ 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 ${\boldsymbol{\beta}}$, not in ${\boldsymbol{x}}$.

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.

1.1. Motivation

If ${\boldsymbol{X}}$ is large and low-rank, one may have many reasons to work with a factorization of ${\boldsymbol{X}}$. We shall discuss three reasons below — algorithmic, infrastructural, and statistical.

Consider data matrices encountered in “recommender systems” in machine learning [2, 14, 20, 26, 27]. 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. ${\boldsymbol{X}}$ 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 ${\boldsymbol{X}}$, 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 [5, 13, 28, 12]. One popular algorithm, alternating-minimization [10], chooses a (small) target rank $k$, and tries to find ${\boldsymbol{U}},{\boldsymbol{V}}$ such that ${\boldsymbol{X}}_{ij}\approx({\boldsymbol{U}}{\boldsymbol{V}})_{ij}$ for all the observed entries $(i,j)$ of ${\boldsymbol{X}}$. As its name suggests, the algorithm alternates between solving for ${\boldsymbol{U}}$ keeping ${\boldsymbol{V}}$ fixed and then solving for ${\boldsymbol{V}}$ keeping ${\boldsymbol{U}}$ fixed. In this case, at no point does the algorithm even form the entire completed (inferred) matrix ${\boldsymbol{X}}$, and the algorithm only has access to factors ${\boldsymbol{U}},{\boldsymbol{V}}$ simply due to algorithmic choices.

There may be other instances where a data scientist may have access to the full matrix ${\boldsymbol{X}}$, but in order to reduce the memory storage footprint, or to communicate the data, may explicitly choose to decompose ${\boldsymbol{X}}\approx{\boldsymbol{U}}{\boldsymbol{V}}$ and discard ${\boldsymbol{X}}$ 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 ${\boldsymbol{X}}$, 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” [15, 29, 18], where one decomposes ${\boldsymbol{X}}$ as the product of two low-rank non-negative matrices ${\boldsymbol{U}},{\boldsymbol{V}}$; the non-negativity is imposed for human interpretability, so that ${\boldsymbol{U}}$ can be interpreted as a documents-by-topics matrix, and ${\boldsymbol{V}}$ as a topics-by-words. In this case, we do not have access to ${\boldsymbol{X}}$ 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 [7]. Hence, the factorization may arise from a purely statistical motivation.

Given a vector ${\boldsymbol{y}}$ (representing age, or document popularity, for example), suppose the data scientist is interested in regressing ${\boldsymbol{X}}$ onto ${\boldsymbol{y}}$, 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 ${\boldsymbol{U}}$ and ${\boldsymbol{V}}$ rather than computing the full system ${\boldsymbol{X}}$?

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.

1.2. Main contribution

We propose two stochastic iterative methods for solving system (1) without computing the product of ${\boldsymbol{U}}$ and ${\boldsymbol{V}}$. 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 (2) and RK iterates to solve (3), that converges to the so-called ordinary least squares solution. We prove linear (“exponential”) convergence to the solution in both cases.

1.3. 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 3 we investigate variations of settings for subsystems (2) and (3) that arise depending on the consistency and size of ${\boldsymbol{X}}$. Section 4 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 5 and conclude with final remarks and future work in Section 6.

1.4. Notation

Here and throughout the paper, matrices and vectors are denoted with boldface letters (uppercase for matrices and lowercase for vectors). We call ${\boldsymbol{X}^{i}}$ the $i^{th}$ row of the matrix ${\boldsymbol{X}}$ and ${\boldsymbol{X}_{(j)}}$ the $j^{th}$ column of ${\boldsymbol{X}}$. The euclidean norm is denoted by $\left\lVert\cdot\right\rVert_{2}$ and the Frobenius norm by $\left\lVert\cdot\right\rVert_{F}$. Lastly, ${\boldsymbol{X}}^{*}$ denotes the adjoint (conjugate transpose) of the matrix ${\boldsymbol{X}}$. Motivated by applications, we allow ${\boldsymbol{X}}$ to be rank deficient and assume that ${\boldsymbol{U}}$ and ${\boldsymbol{V}}$ are full rank.

2. Background and Existing Work

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

2.1. Linear Systems

Linear systems take on one of three settings determined by the size of the system, rank of the matrix ${\boldsymbol{X}}$, and the existence of a solution. First we discuss solutions to systems with full rank matrices ${\boldsymbol{X}}$ 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 (1):

[TABLE]

Clearly, ${\boldsymbol{X}}{\boldsymbol{\beta_{LN}}}={\boldsymbol{y}}$, and all other solutions to an underdetermined system can be written as ${\boldsymbol{b}}={\boldsymbol{\beta_{LN}}}+{\boldsymbol{z}}$ where ${\boldsymbol{X}}{\boldsymbol{z}}=\textbf{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 ${\boldsymbol{X}}$ is full rank, the optimal unique solution is ${\boldsymbol{\beta}}_{\text{uniq}}$ such that ${\boldsymbol{X}}{\boldsymbol{\beta}}_{\text{uniq}}={\boldsymbol{y}}$:

[TABLE]

If there is no exact solution, the system is called an overdetermined inconsistent system. When a system is inconsistent and ${\boldsymbol{X}}$ is full rank, we often seek to minimize the sum of squared residuals, i.e. to find the ordinary least squares solution

[TABLE]

The residual can be written as ${\boldsymbol{r}}={\boldsymbol{X}}{\boldsymbol{\beta_{LS}}}-{\boldsymbol{y}}$. Note that ${\boldsymbol{X}}^{*}{\boldsymbol{r}}=\textbf{0}$, which can be easily seen by substituting ${\boldsymbol{y}}={\boldsymbol{X}}{\boldsymbol{\beta_{LS}}}+{\boldsymbol{r}}$ into (6). For simplicity, we will refer to the matrix ${\boldsymbol{X}}$ of a linear system as consistent or inconsistent when the system itself is consistent or inconsistent.

If the matrix ${\boldsymbol{X}}$ in the linear system ${\boldsymbol{X}}{\boldsymbol{\beta}}={\boldsymbol{y}}$ 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,

[TABLE]

where $(\cdot)^{\dagger}$ is the pseudo-inverse. General solutions to the linear system can be written as ${\boldsymbol{b}}={\boldsymbol{\beta_{LN}}}+{\boldsymbol{z}}$ where ${\boldsymbol{X}}{\boldsymbol{z}}=\textbf{0}$—note that ${\boldsymbol{y}}={\boldsymbol{X}}{\boldsymbol{b}}={\boldsymbol{X}}{\boldsymbol{\beta_{LN}}}+{\boldsymbol{X}}{\boldsymbol{z}}={\boldsymbol{X}}{\boldsymbol{\beta_{LN}}}$. Similar to the full rank case, when the low-rank system is inconsistent, we can write ${\boldsymbol{y}}={\boldsymbol{X}}{\boldsymbol{\beta_{LN}}}+{\boldsymbol{r}}$, again where ${\boldsymbol{X}}^{*}{\boldsymbol{r}}=\textbf{0}$.

2.2. Randomized Kaczmarz and its Extension

The Kaczmarz Algorithm [11] solves a linear system ${\boldsymbol{X}}{\boldsymbol{\beta}}={\boldsymbol{y}}$ by cycling through rows of ${\boldsymbol{X}}$ and projects the estimate onto the solution space given by the chosen row. It was initially proposed by Kaczmarz [11] and has recently regained interest in the setting of computer tomography where it is known as the Algebraic Reconstruction Technique [8, 21, 4, 9]. The randomized variant of the Kaczmarz method introduced by Strohmer and Vershynin [25] was proven to converge linearly in expectation for consistent systems. Formally, given ${\boldsymbol{X}}$ and ${\boldsymbol{y}}$ of (1), RK chooses row $i\in\{1,2,...m\}$ of ${\boldsymbol{X}}$ with probability $\frac{\left\lVert{\boldsymbol{X}^{i}}\right\rVert^{2}_{2}}{\left\lVert{\boldsymbol{X}}\right\rVert^{2}_{F}}$, and projects the previous estimate onto that row with the update

[TABLE]

Needell [22] 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 [30] proposed the Randomized Extended Kaczmarz (REK) algorithm to solve linear systems in all settings. For REK, row $i\in\{1,2,...m\}$ and column $j\in\{1,...n\}$ of ${\boldsymbol{X}}$ are chosen at random with probability

[TABLE]

and starting from ${\boldsymbol{\beta}}_{0}=\textbf{0}$ and ${\boldsymbol{z}}_{0}={\boldsymbol{y}}$, every iteration computes

[TABLE]

REK finds the optimal solution in all linear system settings. In the consistent setting, it behaves as RK. In the overdetermined inconsistent setting, ${\boldsymbol{z}}$ estimates the residual vector ${\boldsymbol{r}}$ and allows ${\boldsymbol{\beta_{t}}}$ 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 [30].

2.3. 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 [3]. Instead of relying on rows of a matrix, the Gauss-Seidel method relies on columns of ${\boldsymbol{X}}$. The randomized variant was studied by Leventhal and Lewis [16] shortly after RK was published. The randomized variant (RGS) requires a column $j$ to be chosen randomly with probability $\frac{\left\lVert{\boldsymbol{X}_{(j)}}\right\rVert^{2}_{2}}{\left\lVert{\boldsymbol{X}}\right\rVert^{2}_{F}}$, and updates at every iteration

[TABLE]

where ${\boldsymbol{e}}_{(j)}$ is the $j^{th}$ basis vector (a vector with 1 in the $j^{th}$ position and 0 elsewhere). Leventhal and Lewis [16] showed that RGS converges linearly in expectation when ${\boldsymbol{X}}$ is overdetermined. However, it fails to find the least norm solution for an underdetermined linear system [19]. 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 ${\boldsymbol{X}}$ exactly as in REK, and then updates at every iteration

[TABLE]

and at any fixed time $t$, outputs ${\boldsymbol{\beta_{t}}}-{\boldsymbol{z}}_{t}$ as the estimated solution to ${\boldsymbol{X}}{\boldsymbol{\beta}}={\boldsymbol{y}}$. Here, ${\boldsymbol{Id}}_{n}$ denotes the $n\times n$ identity matrix. This extension works for all variations of linear systems and was proven to converge linearly in expectation by Ma et al. [19].

The RK and RGS methods along with their extensions are extensively studied and compared in [19]. Table 1 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.

3. Variations of Factored Linear Systems

Our proposed methods rely on interleaving solution estimates to subsystem (2) and subsystem (3). 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 (2) and (3) are determined by ${\boldsymbol{X}}$. In this section, we will discuss when we can expect our methods to solve the full system.

For simplicity in notation, we will denote ${\boldsymbol{\beta_{\star}}}$, $\boldsymbol{x}_{\star}$, and ${\boldsymbol{b}_{\star}}$ as the “optimal” solution of (1), (2), and (3) respectively, as summarized in Table 2. By “optimal” solution for (2) and (3), 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 ${\boldsymbol{U}}{\boldsymbol{V}}$ may be low-rank, ${\boldsymbol{\beta_{\star}}}$ is going to be the least norm solution as described in (7). Table 3 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 3. For this, we split Table 3 into three scenarios : (S1) ${\boldsymbol{U}}$ is overdetermined and consistent, (S2) ${\boldsymbol{U}}$ is underdetermined, (S3a and S3b) ${\boldsymbol{X}}$ 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 ${\boldsymbol{b}_{\star}}={\boldsymbol{\beta_{\star}}}$ for Scenario S1 and ${\boldsymbol{b}_{\star}}\neq{\boldsymbol{\beta_{\star}}}$ 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 (1) affects the consistency of (2) and (3), then show that for Scenario S3a, ${\boldsymbol{b}_{\star}}\neq{\boldsymbol{\beta_{\star}}}$ and for Scenario S3b ${\boldsymbol{b}_{\star}}={\boldsymbol{\beta_{\star}}}$. 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 ${\boldsymbol{U}}$ is overdetermined and consistent, we find that solving (2) and (3) gives us the optimal solution of (1). Indeed, in the case where ${\boldsymbol{V}}$ is overdetermined and consistent, we have:

[TABLE]

In the case where ${\boldsymbol{V}}$ is underdetermined, we have:

[TABLE]

Since ${\boldsymbol{X}}$ is possibly low-rank, we still need to argue that this implies that ${\boldsymbol{b}_{\star}}={\boldsymbol{\beta_{\star}}}$, i.e. ${\boldsymbol{b}_{\star}}$ is indeed the least norm solution to the full system. Suppose towards a contradiction that ${\boldsymbol{\beta_{\star}}}$ is the least norm solution of the full system but not subsystem (3); in other words assume that ${\boldsymbol{\beta_{\star}}}={\boldsymbol{b}_{\star}}+{\boldsymbol{b}}$ where ${\boldsymbol{V}}{\boldsymbol{b}}=\textbf{0}$ and ${\boldsymbol{b}}$ is nontrivial. Multiplying both sides by ${\boldsymbol{X}}$, we see that since ${\boldsymbol{\beta_{\star}}}$ has a nontrivial component ${\boldsymbol{b}}$ such that ${\boldsymbol{X}}{\boldsymbol{b}}=\textbf{0}$, it cannot be the least norm solution to the full system as assumed, reaching a contradiction. Therefore, in the consistent case when ${\boldsymbol{U}}$ is overdetermined, we have that ${\boldsymbol{b}_{\star}}={\boldsymbol{\beta_{\star}}}$, and may hope that our proposed methods will be able to solve the full system (1) utilizing the subsystems (2) and (3).

• •

Scenario S2: U underdetermined. When ${\boldsymbol{U}}$ is underdetermined, solving (2) and (3) for their optimal solutions does not guarantee the optimal solution of the full system. Intuitively, (2) has infinitely many solutions and $\boldsymbol{x}_{\star}={\boldsymbol{x}}_{LN}\not={\boldsymbol{V}}{\boldsymbol{\beta_{\star}}}$. Mathematically, investigating ${\boldsymbol{b}_{\star}}$, we find that

[TABLE]

where we rewrite $\boldsymbol{x}_{\star}={\boldsymbol{V}}{\boldsymbol{b}_{\star}}+{\boldsymbol{r}}_{V}$ for ${\boldsymbol{r}}_{V}\in null({\boldsymbol{V}}^{*})$ since subsystem (3) may be inconsistent. Note that if subsystem (3) is consistent then we simply have ${\boldsymbol{r}}_{V}=0$ and the above calculation still carries through. Therefore, we do not expect our proposed methods to succeed when ${\boldsymbol{U}}$ 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 ${\boldsymbol{X}}$ being inconsistent implies about the subsystems (2) and (3). In particular, one needs to determine whether inconsistency in the full system creates inconsistencies in the individual subsystems. If ${\boldsymbol{X}}$ is inconsistent then we have ${\boldsymbol{X}}{\boldsymbol{\beta_{\star}}}+{\boldsymbol{r}}={\boldsymbol{y}}$ where ${\boldsymbol{\beta_{\star}}}$ is the optimal solution of (1) and ${\boldsymbol{X}}^{*}{\boldsymbol{r}}=\textbf{0}$. Now, consider decomposing ${\boldsymbol{r}}={\boldsymbol{r}}_{1}+{\boldsymbol{r}}_{2}$ where ${\boldsymbol{U}}^{*}{\boldsymbol{r}}_{1}=\textbf{0}$, ${\boldsymbol{U}}^{*}{\boldsymbol{r}}_{2}\neq\textbf{0}$, and ${\boldsymbol{V}}^{*}{\boldsymbol{U}}^{*}{\boldsymbol{r}}_{2}=0$. Notice that ${\boldsymbol{X}}^{*}{\boldsymbol{r}}={\boldsymbol{V}}^{*}{\boldsymbol{U}}^{*}({\boldsymbol{r}}_{1}+{\boldsymbol{r}}_{2})={\boldsymbol{V}}^{*}{\boldsymbol{U}}^{*}{\boldsymbol{r}}_{1}+{\boldsymbol{V}}^{*}{\boldsymbol{U}}^{*}{\boldsymbol{r}}_{2}=\textbf{0}$, as desired. We want to decompose the full system ${\boldsymbol{X}}{\boldsymbol{\beta_{\star}}}+{\boldsymbol{r}}={\boldsymbol{y}}$ into two subsystems. Following a similar thought process as before, we choose to decompose our full system into the following:

[TABLE]

Clearly, (9) is inconsistent since ${\boldsymbol{U}}^{*}{\boldsymbol{r}}_{1}=\textbf{0}$ and ${\boldsymbol{U}}^{*}{\boldsymbol{r}}_{2}\neq\textbf{0}$.

Because ${\boldsymbol{U}}$ must be overdetermined for ${\boldsymbol{X}}$ to be inconsistent, $\boldsymbol{x}_{\star}=({\boldsymbol{U}}^{*}{\boldsymbol{U}})^{-1}{\boldsymbol{U}}^{*}({\boldsymbol{y}}-{\boldsymbol{r}}_{1}-{\boldsymbol{r}}_{2})=({\boldsymbol{U}}^{*}{\boldsymbol{U}})^{-1}{\boldsymbol{U}}^{*}({\boldsymbol{y}}-{\boldsymbol{r}}_{2})$ is the least squares solution to (9).

• –

Case S3a: ${\boldsymbol{V}}$ overdetermined. Note that the second subsystem (10) is possibly inconsistent (since there may be a component $\boldsymbol{x}_{\star}$ of in the null space of ${\boldsymbol{V}}^{*}$). Writing $\boldsymbol{x}_{\star}={\boldsymbol{V}}{\boldsymbol{b}_{\star}}+{\boldsymbol{r}}_{V}$ such that ${\boldsymbol{r}}_{V}\in null({\boldsymbol{V}}^{*})$, we have

[TABLE]

Similar to Scenario S2, if subsystem (10) is indeed consistent then ${\boldsymbol{r}}_{V}=0$ and the above calculation still carries through. Therefore, in this case, we do not expect to find the optimal solution to (1).

• –

Case S3b: ${\boldsymbol{V}}$ underdetermined. In this case, ${\boldsymbol{r}}_{2}=\textbf{0}$ and solving (9) and (10) obtains the optimal solution to the full system since

[TABLE]

Following the same argument as in Scenario S1 when ${\boldsymbol{V}}$ is underdetermined, we reach the conclusion that ${\boldsymbol{b}_{\star}}={\boldsymbol{\beta_{\star}}}$. *Thus, in this case our methods have the potential to solve the full system.

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

4. Methods and Main Results

Our approach intertwines two iterative methods to solve subsystem (2) followed by subsystem (3). For the consistent setting, we propose Algorithm 1 which uses an iterate of RK on (2) intertwined with an iterate of RK to solve (3). For the inconsistent setting, we propose using REK to solve subsystem (9) followed by RK to solve subsystem (10) as shown in Algorithm 2. We view the latter method as a more practical approach, and the former as interesting from a theoretical point of view. Recall that $\boldsymbol{x}_{t}^{p}$ is the $p^{th}$ element in the vector $\boldsymbol{x}_{t}$. 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 (2) and (3); this has a couple of advantages over solving each subsystem separately. First, if we are given some tolerance $\epsilon$ 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.

4.1. Main result

Our main result shows that Algorithm 1 and Algorithm 2 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 ${\boldsymbol{A}}$ we write

[TABLE]

where $\sigma_{\min}^{2}({\boldsymbol{A}})$ is the smallest non-zero singular value of ${\boldsymbol{A}}$, and ${\kappa_{A}^{2}}$ is the squared condition number of ${\boldsymbol{A}}$. 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.

Theorem 1**.**

Let ${\boldsymbol{X}}$ be low rank, ${\boldsymbol{X}}={\boldsymbol{U}}{\boldsymbol{V}}$ such that ${\boldsymbol{U}}\in\mathbb{C}^{m\times k}$ and ${\boldsymbol{V}}\in\mathbb{C}^{k\times n}$ are full rank, and the systems ${\boldsymbol{X}}{\boldsymbol{\beta}}={\boldsymbol{y}}$, ${\boldsymbol{U}}{\boldsymbol{x}}={\boldsymbol{y}}$, and ${\boldsymbol{V}}{\boldsymbol{b}}={\boldsymbol{x}}$ have optimal solutions ${\boldsymbol{\beta_{\star}}}$, $\boldsymbol{x}_{\star}$, and ${\boldsymbol{b}_{\star}}$ respectively, and $\alpha_{U}$, $\alpha_{V},\kappa^{2}_{U}$ are as defined in (11) and (12). Setting ${\boldsymbol{b}}_{0}=\textbf{0}$ and assuming $k<m,n$, we have

1. (a)

if ${\boldsymbol{X}}{\boldsymbol{\beta}}={\boldsymbol{y}}$ is consistent, then ${\boldsymbol{b}_{\star}}={\boldsymbol{\beta_{\star}}}$ and Algorithm 1 converges with expected error

[TABLE]

where $\alpha_{\text{max}}=\max\{\alpha_{U},\alpha_{V}\}$ and $\gamma_{1}=\min\{\frac{\alpha_{U}}{\alpha_{V}},\frac{\alpha_{V}}{\alpha_{U}}\}$. 2. (b)

if ${\boldsymbol{X}}{\boldsymbol{\beta}}={\boldsymbol{y}}$ is inconsistent, then ${\boldsymbol{b}_{\star}}={\boldsymbol{\beta_{\star}}}$ and Algorithm 2 converges with expected error

[TABLE]

where $\tilde{\alpha}_{\text{max}}=\max\{\sqrt{\alpha_{U}},\alpha_{V}\}$ and $\gamma_{2}=\min\{\frac{\sqrt{\alpha_{U}}}{\alpha_{V}},\frac{\alpha_{V}}{\sqrt{\alpha_{U}}}\}$.

Remarks.

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

2. Empirically, our experiments in the next section suggest that ${\boldsymbol{U}}$ and ${\boldsymbol{V}}$ can be substantially better conditioned than ${\boldsymbol{X}}$.

3. Algorithm 1 is interesting to discuss from a theoretical standpoint but in applications Algorithm 2 is more practical as linear systems are typically inconsistent. Algorithm 1 can be utilized in applications if error in the solution is tolerable. In particular, if ${\boldsymbol{X}}{\boldsymbol{\beta}}={\boldsymbol{y}}$ is inconsistent then Algorithm 1 will converge in expectation to some convergence horizon. This can be see by replacing the use of Proposition 1 in the bound (14) with the convergence bound of RK on inconsistent linear systems found in Theorem 2.1 of [22].

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 ${\boldsymbol{U}}$ and ${\boldsymbol{V}}$, 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.

4.2. 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 ${\mathbb{E}}_{b}$ denote the expected value taken over the choice of rows in ${\boldsymbol{V}}$ and ${\mathbb{E}}_{x}$ the expected value taken over the choice of rows in ${\boldsymbol{U}}$ and when necessary the choice of columns in ${\boldsymbol{U}}$. Also, let ${\mathbb{E}}$ denote the full expected value (over all random variables and iterations) and ${\mathbb{E}}^{t-1}$ be the expectation conditional on the first $t-1$ iterations.

Proposition 1**.**

([24, Theorem 2]) Given a consistent linear system ${\boldsymbol{X}}{\boldsymbol{\beta}}={\boldsymbol{y}}$, the Randomized Kaczmarz algorithm, with initialization ${\boldsymbol{\beta}}_{0}=\textbf{0}$, as described in Section 2.2 converges to the optimal solution ${\boldsymbol{\beta_{\star}}}$ with expected error

[TABLE]

where $\alpha_{X}$ is as defined in (11).

Proposition 2**.**

([30, Theorem 8]) Given a linear system ${\boldsymbol{X}}{\boldsymbol{\beta}}={\boldsymbol{y}}$, the Randomized Extended Kaczmarz algorithm, with initialization ${\boldsymbol{\beta}}_{0}=\textbf{0}$, as described in Section 2.2 converges to the optimal solution ${\boldsymbol{\beta_{\star}}}$ with expected error

[TABLE]

where $\alpha_{X}$ is as defined in (11) and ${\kappa_{X}^{2}}$ is as defined in (12).

The proof of Theorem 1 builds directly on two useful lemmas. Lemma 1 addresses the impact of intertwining the algorithms. In particular, it shows useful relationships involving ${\boldsymbol{\tilde{b}}_{t}}$, the RK update solving the linear system ${\boldsymbol{V}}{\boldsymbol{b}}=\boldsymbol{x}_{\star}$ at the $t^{th}$ iteration (with ${\boldsymbol{b}_{t-1}}$ as the previous estimate), and our update ${\boldsymbol{b}_{t}}$. Lemma 2 states that conditional on the first $t-1$ iterations, we can split the norm squared error $\|{\boldsymbol{b}_{t}}-{\boldsymbol{b}_{\star}}\|^{2}$ into two terms relating to the error from solving subsystem (2) and the error from solving subsystem (3). To complete the proof of Theorem 1, we bound the error from solving (2) depending on whether we use RK (as in Algorithm 1) or REK (as in Algorithm 2) then apply the law of iterated expectations to bound the error from solving (3). We now state the aforementioned lemmas, and then formally prove the theorem.

Lemma 1**.**

Let ${\boldsymbol{\tilde{b}}_{t}}={\boldsymbol{b}_{t-1}}+\frac{\boldsymbol{x}_{\star}^{p}-{\boldsymbol{V}^{p}}{\boldsymbol{b}_{t-1}}}{\|{\boldsymbol{V}^{p}}\|^{2}}({\boldsymbol{V}^{p}})^{*}$. In Algorithm 1 and Algorithm 2 we have that:

1. (a)

${\mathbb{E}}_{b}^{t-1}\left\langle{\boldsymbol{b}_{t}}-{\boldsymbol{\tilde{b}}_{t}},{\boldsymbol{\tilde{b}}_{t}}-{\boldsymbol{b}_{\star}}\right\rangle=0,$** 2. (b)

$\|{\boldsymbol{\tilde{b}}_{t}}-{\boldsymbol{b}_{\star}}\|^{2}=\|{\boldsymbol{b}_{t-1}}-{\boldsymbol{b}_{\star}}\|^{2}-\|{\boldsymbol{\tilde{b}}_{t}}-{\boldsymbol{b}_{t-1}}\|^{2}.$**

In words, part (a) states that the difference between an RK iterate solving the exact linear system ${\boldsymbol{V}}{\boldsymbol{b}}=\boldsymbol{x}_{\star}$ and our RK iterate (which solves the linear system resulting from intertwining ${\boldsymbol{V}}{\boldsymbol{b}}=\boldsymbol{x}_{t}$), is orthogonal to ${\boldsymbol{\tilde{b}}_{t}}-{\boldsymbol{b}_{\star}}$. This will come in handy in Lemma 2. Part (b) is a Pythagoras-style statement, which follows from well-known orthogonality properties of RK updates, included here for simplicity and completeness.

Proof.

To prove statement (b), we note that $({\boldsymbol{\tilde{b}}_{t}}-{\boldsymbol{b}_{t-1}})$ is parallel to ${\boldsymbol{V}^{p}}$ and $({\boldsymbol{\tilde{b}}_{t}}-{\boldsymbol{b}_{\star}})$ is perpendicular to ${\boldsymbol{V}^{p}}$ since ${\boldsymbol{V}^{p}}({\boldsymbol{\tilde{b}}_{t}}-{\boldsymbol{b}_{\star}})={\boldsymbol{V}^{p}}({\boldsymbol{b}_{t-1}}+\frac{\boldsymbol{x}_{\star}^{p}-{\boldsymbol{V}^{p}}{\boldsymbol{b}_{t-1}}}{\|{\boldsymbol{V}^{p}}\|^{2}}({\boldsymbol{V}^{p}})^{*}-{\boldsymbol{b}_{\star}})={\boldsymbol{V}^{p}}{\boldsymbol{b}_{t-1}}+\boldsymbol{x}_{\star}^{p}-{\boldsymbol{V}^{p}}{\boldsymbol{b}_{t-1}}-\boldsymbol{x}_{\star}^{p}=\textbf{0}$. We apply the Pythagorean Theorem to obtain the desired result.

We prove statement (a) by direct substitution and expansion, as follows:

[TABLE]

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 ${\boldsymbol{V}}$. In step ($iv$), we use the fact that for $k<m,n$, subsystem (3) is always consistent (since ${\boldsymbol{V}}$ is underdetermined) to make the substitution $\boldsymbol{x}_{\star}={\boldsymbol{V}}{\boldsymbol{b}_{\star}}$.

∎

Lemma 2**.**

In Algorithm 1 and Algorithm 2, we can bound the expected norm squared error of ${\boldsymbol{b}_{t}}-{\boldsymbol{b}_{\star}}$ as

[TABLE]

We investigate the expectation of the norm squared error of ${\boldsymbol{b}_{t}}-{\boldsymbol{b}_{\star}}$ conditional on the first $t-1$ iterations and over the choice of rows of ${\boldsymbol{V}}$. We keep ${\mathbb{E}}_{x}^{t-1}$ in our bound as this expectation will depend on whether Algorithm 1 or Algorithm 2 is being used.

Proof.

[TABLE]

Steps ($iii$) and ($iv$) are applications of Lemma 1(a) and Lemma 1(b) respectively, and step ($v$) follows from the definition of each term and simplification using the fact that ${\boldsymbol{V}}{\boldsymbol{b}_{\star}}=\boldsymbol{x}_{\star}$.

Now, we evaluate the conditional expectation on the choices of rows of ${\boldsymbol{V}}$ to complete the proof:

[TABLE]

In step ($vi$), we use iterated expectations to split the expected value ${\mathbb{E}}^{t-1}={\mathbb{E}}_{x}^{t-1}{\mathbb{E}}_{b}^{t-1}$. Step $(vii)$ uses the fact that $\|{\boldsymbol{V}}({\boldsymbol{b}_{t-1}}-{\boldsymbol{b}_{\star}})\|^{2}\geq\sigma_{\min}^{2}({\boldsymbol{V}})\|{\boldsymbol{b}_{t-1}}-{\boldsymbol{b}_{\star}}\|^{2}$ since ${\boldsymbol{b}_{t-1}}-{\boldsymbol{b}_{\star}}$ are in the row span of ${\boldsymbol{V}}$ for all $t$. We simplify and obtain the desired bound. ∎

4.3. Proof of main result

We now have all the ingredients we need to prove Theorem 1, which we now proceed to below.

Proof of Theorem 1.

The fact that ${\boldsymbol{b}_{\star}}={\boldsymbol{\beta_{\star}}}$ 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 1, we only need to invoke the statement of Lemma 2 and bound the term ${\mathbb{E}}_{x}^{t-1}\|\boldsymbol{x}_{t}-\boldsymbol{x}_{\star}\|^{2}$ using Proposition 1 or Proposition 2 depending on whether we are using Algorithm 1 or Algorithm 2, respectively.

1. (a)

For Algorithm 1, plugging Proposition 1 into the statement of Lemma 2 yields

[TABLE]

Let $\alpha_{\text{max}}=\max\{\alpha_{U},\alpha_{V}\}$, $\gamma_{1}=\min\{\frac{\alpha_{U}}{\alpha_{V}},\frac{\alpha_{V}}{\alpha_{U}}\}$, and note that $\gamma_{1}<1$ if $\alpha_{U}\neq\alpha_{V}$. Taking expectations over the randomness from the first $t-1$ iterations and using the Law of Iterated Expectation, we have

[TABLE]

where step (i) uses the fact that the summation is symmetric with respect to $\alpha_{U}$ and $\alpha_{V}$. Step (ii) uses the fact that $\alpha_{\text{max}}^{t-1}\alpha_{U}=\alpha_{\text{max}}^{t}$ if $\alpha_{\text{max}}=\alpha_{U}$ and $\alpha_{\text{max}}^{t-1}\alpha_{U}<\alpha_{\text{max}}^{t}$ if $\alpha_{\text{max}}=\alpha_{V}$. Now, if $\alpha_{U}=\alpha_{V}=\alpha_{\text{max}}$ then we have

[TABLE]

where the second term of $(ii)$ approaches 0 as $t\rightarrow\infty$ since $\alpha_{\text{max}}<1$. 2. (b)

For Algorithm 2, plugging Proposition 2 into the statement of Lemma 2 yields

[TABLE]

Taking expectations over the remaining randomness, and using the fact that $\alpha_{U}^{\lfloor\frac{t-h}{2}\rfloor}\leq\alpha_{U}^{\frac{t-h-1}{2}}$ since $\alpha_{U}<1$, we have

[TABLE]

Using the same techniques as in (a), let $\tilde{\alpha}_{\text{max}}=\max\{\alpha_{V},\sqrt{\alpha_{U}}\}$ and $\gamma_{2}=\min\{\frac{\alpha_{V}}{\sqrt{\alpha_{U}}},\frac{\sqrt{\alpha_{U}}}{\alpha_{V}}\}$ and again noting $\gamma_{2}<1$ when $\alpha_{V}\neq\sqrt{\alpha_{U}}$, we can write:

[TABLE]

When $\alpha_{V}=\sqrt{\alpha_{U}}=\tilde{\alpha}_{\text{max}}$,

[TABLE]

where the second term of $(i)$ goes to 0 as $t$ goes to infinity.

This concludes the proof of the theorem. ∎

5. 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 (2) and ALG2 is the algorithm used to solve subsystem (3). When an algorithm’s name is used alone, we imply applying the algorithm on the full system (1).

In Figure 1 we show our first set of experiments. Entries of ${\boldsymbol{U}}$, ${\boldsymbol{V}}$, and ${\boldsymbol{\beta}}$ are drawn from a standard Gaussian distribution. We set ${\boldsymbol{X}}={\boldsymbol{U}}{\boldsymbol{V}}$ and ${\boldsymbol{y}}={\boldsymbol{X}}{\boldsymbol{\beta}}$ if ${\boldsymbol{X}}$ is consistent and ${\boldsymbol{y}}={\boldsymbol{X}}{\boldsymbol{\beta}}+{\boldsymbol{r}}$ where ${\boldsymbol{r}}\in null({\boldsymbol{X}}^{*})$ (computed in Matlab using null() function) if ${\boldsymbol{X}}$ is inconsistent. In this first set of experiments, $m,n,k\in\{100,\,150,\,200\}$ depending on the desired size of $k$ with respect to the over or underdetermined-ness of ${\boldsymbol{X}}$. For example, if $k<m,n$ and ${\boldsymbol{X}}$ is overdetermined then $k=100$, $m=200$, and $n=150$. The plots show iteration vs $\ell_{2}$-error, $\|{\boldsymbol{b}_{t}}-{\boldsymbol{\beta_{\star}}}\|^{2}$, of each method averaged over 40 runs and allowing each algorithm to run $7\times 10^{4}$ iterations. The layout of Figure 1 is exactly as in Table 3. For each row, we have a different setting for ${\boldsymbol{X}}$ and for each column, we vary the size of $k$ depending on the size of ${\boldsymbol{X}}$. Looking at the overall trends, we see that when $k<m,n$ and when ${\boldsymbol{X}}$ 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 3. Looking at each individual subplot, we also find what one would expect according to Table 1. In other words, if ${\boldsymbol{U}}$ or ${\boldsymbol{V}}$ is in one of the settings where RK or RGS are expected to fail then RK-RK or RGS-RGS fail as well.

When ${\boldsymbol{X}}$ is overdetermined, inconsistent and $k<m,n$ we have that ${\boldsymbol{V}}$ 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 2 has the same set up as discussed in the previous experiment with the exception of using larger random matrices with ${\boldsymbol{X}}:1200\times 750$ and $k=500$. In Figure 2(a) we plot iteration vs $\ell_{2}$-error and in Figure 2(b) we plot FLOPS vs $\ell_{2}$-error. The errors are averaged over 40 runs, with shaded regions representing error within two standard deviations from the mean. Note that Algorithm 2 performs excellently in practice, better than our theoretical upper bound as computed in Theorem 1. 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 $\ell_{2}$-error at the $t^{th}$ iteration over 40 runs and shaded regions representing the $\ell_{2}$-error within two standard deviations. In addition to empirical performance, we also plot the theoretical convergence bound derived in Theorem 1 (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 [17]. 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 ${\boldsymbol{U}}$ and ${\boldsymbol{V}}$ using Matlab’s nnmf() function for nonnegative matrix factorization (recall the motivations from the first section). Figure 3(a) shows the results from this experiment. The conditioning of ${\boldsymbol{X}}$, ${\boldsymbol{U}}$, and ${\boldsymbol{V}}$ are ${\kappa_{X}^{2}}=2.46\times 10^{3}$, ${\kappa_{U}^{2}}=25.96$, and ${\kappa_{V}^{2}}=4.20$ respectively. We plot the $\ell_{2}$-error averaged over 40 runs. Since ${\boldsymbol{X}}$ has such a large condition number, this impacts the convergence of REK on ${\boldsymbol{X}}$ 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 ${\boldsymbol{U}}$ and ${\boldsymbol{V}}$.

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 ${\boldsymbol{U}}$ and ${\boldsymbol{V}}$ in the same way as with the wine data set. Figure 3(b) shows the results from this experiment. The conditioning of ${\boldsymbol{X}}$, ${\boldsymbol{U}}$, and ${\boldsymbol{V}}$ are ${\kappa_{X}^{2}}=94.27$, ${\kappa_{U}^{2}}=54.91$, and ${\kappa_{V}^{2}}=2.99$ respectively. Similar to Figure 3(a), we see that the convergence of REK suffers from the poorly conditioned matrix ${\boldsymbol{X}}$. 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 ${\boldsymbol{U}}:10^{6}\times 10^{3}$ and ${\boldsymbol{V}}:10^{3}\times 10^{4}$. These matrices are so large that the matrix product ${\boldsymbol{U}}{\boldsymbol{V}}$ cannot be computed in Matlab due to memory constraints. These results are shown in Figure 4. We see that without needing to do the matrix computation, we are still able to find the solution to the linear system ${\boldsymbol{U}}{\boldsymbol{V}}{\boldsymbol{\beta}}={\boldsymbol{y}}$.

Lastly, we present the performance of our methods on a large real world data set. We use the Yelp challenge data set [1]. In our setting, we let ${\boldsymbol{X}}:10^{5}\times 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 ${\boldsymbol{X}}$ 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, ${\boldsymbol{y}}$ is a vector that represents the number of stars a review received. We choose $k=5000$. Figure 4(b) shows the results from this experiment using REK, REK-REK, and REK-RK. The conditioning of ${\boldsymbol{X}}$, ${\boldsymbol{U}}$, and ${\boldsymbol{V}}$ are ${\kappa_{X}^{2}}=127.3592$, ${\kappa_{U}^{2}}=24.274$, and ${\kappa_{V}^{2}}=19.096$ respectively. In this large real world data set, we can again see the usefulness of our proposed methods when we are given ${\boldsymbol{X}}={\boldsymbol{U}}{\boldsymbol{V}}$.

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.

6. 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 3. 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.