Total Least Squares Regression in Input Sparsity Time
Huaian Diao Zhao Song David P. Woodruff Xin Yang









  The Infall Motion in the Low-Mass Protostellar Binary NGC 1333 IRAS 4A1/4A2
    I-Ta Hsieh
     
===========================================================================


In the total least squares problem, one is given an m × n matrix A, and an m × d matrix B, and one seeks to “correct” both A and B, obtaining matrices Â and B̂, so that there exists an X satisfying the equation ÂX = B̂. Typically the problem is overconstrained, meaning that m ≫max(n,d). The cost of the solution Â, B̂ is given by A-Â_F^2 + B - B̂_F^2. We give an algorithm for finding a solution X to the linear system ÂX=B̂ for which the cost A-Â_F^2 + B-B̂_F^2 is at most a multiplicative (1+ϵ) factor times the optimal cost, up to an additive error η that may be an arbitrarily small function of n. Importantly, our running time is Õ( (A) + (B) ) + (n/ϵ) · d, where for a matrix C, (C) denotes its number of non-zero entries. Importantly, our running time does not directly depend on the large parameter m. As total least squares regression is known to be solvable via low rank approximation, a natural approach is to invoke fast algorithms for approximate low rank approximation, obtaining matrices Â and B̂ from this low rank approximation, and then solving for X so that ÂX = B̂. However, existing algorithms do not apply since in total least squares the rank of the low rank approximation needs to be n, and so the running time of known methods would be at least mn^2. In contrast, we are able to achieve a much faster running time for finding X by never explicitly forming the equation Â X = B̂, but instead solving for an X which is a solution to an implicit such equation.
Finally, we generalize our algorithm to the total least squares problem with regularization. 




  
  empty











§ INTRODUCTION







In the least squares regression problem, we are given an m × n matrix A and an m × 1 vector b, and we seek to find an x ∈ℝ^n which minimizes Ax-b_2^2. A natural geometric interpretation is that there is an unknown hyperplane in ℝ^n+1, specified by the normal vector x, for which we have m points on this hyperplane, the i-th of which is given by (A_i, ⟨ A_i, x ⟩), where A_i is the i-th row of A. However, due to noisy observations, we do not see the points (A_i, ⟨ A_i, x ⟩), but rather only see the point (A_i, b_i), and we seek to find the hyperplane which best fits these points, where we measure (squared) distance only on the (n+1)-st coordinate. This naturally generalizes to the setting in which B is an m × d matrix, and the true points have the form (A_i, A_i X) for some unknown n × d matrix X. This setting is called multiple-response regression, in which one seeks to find X to minimize AX-B_F^2, where for a matrix Y, Y_F^2 is its squared Frobenius norm, i.e., the sum of squares of each of its entries. This geometrically corresponds to the setting when the points live in a lower n-dimensional flat of ℝ^n+d, rather than in a hyperplane. 

While extremely useful, in some settings the above regression model may not be entirely realistic. For example, it is quite natural that the matrix A may also have been corrupted by measurement noise. In this case, one should also be allowed to first change entries of A, obtaining a new m × n matrix Â, then try to fit B to Â by solving a multiple-response regression problem. One should again be penalized for how much one changes the entries of A, and this leads to a popular formulation known as the total least squares optimization problem min_Â, XA-Â_F^2 + ÂX-B_F^2. Letting C = [A, B], one can more compactly write this objective as min_Ĉ = [Â, B̂]C-Ĉ_F^2, where it is required that the columns of B̂ are in the column span of Â. Total least squares can naturally capture many scenarios that least squares cannot. For example, imagine a column of B is a large multiple λ· a of a column a of A that has been corrupted and sent to 0. Then in least squares, one needs to pay λ^2 a_2^2, but in total least squares one can “repair A" to contain the column a, and just pay a_2^2. We refer the reader to <cit.> for an overview of total least squares. There is also
a large amount of work on total least squares with regularization <cit.>. 

Notice that Ĉ has rank n, and therefore the optimal cost is at least C-C_n_F^2, where C_n is the best rank-n approximation to C. If, in the optimal rank-n approximation C_n, one has the property that the last d columns are in the column span of the first n columns, then the optimal solution Ĉ to total least squares problem is equal to C_n, and so the total least squares cost is the cost of the best rank-n approximation to C. In this case, and only in this case, there is a closed-form solution. However, in general, this need not be the case, and C-C_n_F^2 may be strictly smaller than the total least squares cost. Fortunately, though, it cannot be much smaller, since one can take the first n columns of C_n, and for each column that is not linearly independent of the remaining columns, we can replace it with an arbitrarily small multiple of one of the last d columns of C_n which is not in the span of the first n columns of C_n. Iterating this procedure, we find that there is a solution to the total least squares problem which has cost which is arbitrarily close to C-C_n_F^2. We describe this procedure in more detail below.

The above procedure of converting a best rank-n approximation to an arbitrarily close solution to the total least squares problem can be done efficiently given C_n, so this shows one can get an arbitrarily good approximation by computing a truncated singular value decompostion (SVD), which is a standard way of solving for C_n in O(m(n+d)^2) time. However, given the explosion of large-scale datasets these days, this running time is often prohibitive, even for the simpler problem of multiple response least squares regression. Motivated by this, an emerging body of literature has looked at the sketch-and-solve paradigm, where one settles for randomized approximation algorithms which run in much faster, often input sparsity time. Here by input-sparsity, we mean in time linear in the number (C) of non-zero entries of the input description C = [A, B]. By now, it is known, for example, how to output a solution matrix X to multiple response least squares regression satisfying AX-B_F^2 ≤ (1+ϵ) min_X'AX'-B_F^2, in (A) + (B)+ (nd/ϵ) time. This algorithm works for arbitrary input matrices A and B, and succeeds with high probability over the algorithm's random coin tosses. For a survey of this and related results, we refer the reader to <cit.>. 

Given the above characterization of total least squares as a low rank approximation problem, it is natural to ask if one can directly apply sketch-and-solve techniques to solve it. Indeed, for low rank approximation, it is known how to find a rank-k matrix Ĉ for which C-Ĉ_F^2 ≤ (1+ϵ)C-C_k_F^2 in time (C) + m · k^2/ϵ <cit.>, using the fastest known results. Here, recall, we assume m ≥ n+d. From an approximation point of view, this is fine for the total least squares problem, since this means after applying the procedure above to ensure the last d columns of Ĉ are in the span of the first n columns, and setting k = n in the low rank approximation problem, our cost will be at most (1+ϵ)C-C_n_F^2 + η, where η can be made an arbitrarily small function of n. Moreover, the optimal total least squares cost is at least C-C_n_F^2, so our cost is a (1+ϵ)-relative error approximation, up to an arbitarily small additive η. 

Unfortunately, this approach is insufficient for total least squares, because in the total least squares problem one sets k = n, and so the running time for approximate low rank approximation becomes (C) + m · n^2/ϵ. Since n need not be that small, the m · n^2 term is potentially prohibitively large. Indeed, if d ≤ n, this may be much larger than the description of the input, which requires at most mn parameters. Note that just outputting
Ĉ may take m · (n+d) parameters to describe. However, as in
the case of regression, one is often just interested in the matrix X or
the hyperplane x for ordinary least squares regression. Here
the matrix X for total least squares can be described using only nd
parameters, and so one could hope for a much faster running time.



 §.§ Our Contributions


Our main contribution is to develop a (1+ϵ)-approximation to the total least squares regression problem, returning a matrix X ∈ℝ^n × d for which there exist Â∈ℝ^m × n and
B̂∈ℝ^m × d for which Â X = B̂ and
C-Ĉ_F^2 ≤ (1+ϵ)C-C_n_F^2 + η, where
C = [A, B], Ĉ = [Â, B̂], C_n is the best rank-n
approximation to C, and η is an arbitrarily small function of n. 
Importantly, we achieve a running time of 
Õ((A) + (B)) + (n/ϵ) · d.

Notice that this running time may be faster than the time it takes even to write down Â and B̂. 
Indeed, although one can write A and B down in (A) + (B) time, it could be that 
the algorithm can only efficiently find an Â and a B̂ that are dense; nevertheless the algorithm does
not need to write such matrices down, as it is only interested in outputting the solution X to the equation
ÂX = B̂. This is motivated by applications in which one wants generalization error. Given X, and a 
future y ∈ℝ^n, one can compute yX to predict the remaining unknown d coordinates of the extension
of y to n+d dimensions. 

Our algorithm is inspired by using dimensionality reduction techniques for low rank approximation, such as fast
oblivious “sketching” matrices, as well as leverage score sampling. The rough idea is to
quickly reduce the low rank approximation problem to a problem of the form
min_rank - n  Z ∈ℝ^d_1 × s_1(D_2 C D_1) Z (S_1 C) - D_2 C_F, where d_1, s_1 = O(n/ϵ), D_2
and D_1 are row and column subset selection matrices, and S_1 is a so-called CountSketch matrix, which is a fast
oblivious projection matrix. We describe the matrices D_1, D_2, and S_1 in more detail in the next section, though
the key takeaway message is that (D_2 C D_1), (S_1 C), and (D_2 C) are each efficiently computable small matrices 
with a number of non-zero entries no larger than that of C. Now the problem is a small, rank-constrained regression
problem for which there are closed form solutions for Z. We then need additional technical work, of the form described above, in order to 
find an X ∈ℝ^n × d given Z, and to ensure that X is the solution to an equation of the form
Â X = B̂. 
Surprisingly, fast sketching methods have not been applied 
to the total least squares
problem before, and we consider this application to be one of the main 
contributions of this paper. 

We carefully bound the running time at each step to achieve 
Õ((A) + (B) + (n/ϵ) d) 
overall time, and prove its overall approximation ratio. 
Our main result is Theorem <ref>.










We also generalize the theorem to the important case of total least squares regression with regularization; see Theorem <ref> for a precise statement. 

We empirically validate our algorithm on real and synthetic data sets. As expected, on a number of datasets 
the total least squares error can be much smaller than the error of ordinary least squares regression. We then implement
our fast total least squares algorithm, and show it is roughly 20-40 times faster than computing the exact solution to
total least squares, while retaining 95% accuracy. 




Notation. 
For a function f, we define O(f) to be f·log^O(1)(f). 
For vectors x,y∈^n, let ⟨ x,y⟩:=∑_i=1^n x_iy_i denote the inner product of x and y. Let (A) denote the number of nonzero entries of A. Let (A) denote the determinant of a square matrix A. Let A^⊤ denote the transpose of A. Let A^† denote the Moore-Penrose pseudoinverse of A. Let A^-1 denote the inverse of a full rank square matrix. Let A_F denote the Frobenius norm of a matrix A, i.e., A_F= (∑_i ∑_j A_i,j^2)^1/2.


Sketching matrices play an important role in our algorithm. 
Their usefulness will be further explained in Section <ref>.
The reader can refer to Appendix <ref> for detailed introduction.






§ PROBLEM FORMULATION


We first give the precise definition of the exact (i.e., non-approximate) version of the total least squares problem, and then define the approximate case.

Given two matrices A ∈^m × n, B ∈^m × d, let C = [A,   B] ∈^m × ( n + d ). The goal is to solve the following minimization problem:

    min_X ∈^n × d, Δ A ∈^m × n, Δ B ∈^m × d     [ Δ A,   Δ B ] _F
     subject to      ( A + Δ A ) X = ( B + Δ B )


It is known that total least squares problem has a closed form solution. 
For a detailed discussion, see Appendix <ref>.
It is natural to consider the approximate version of total least squares:

Given two matrices A ∈^m × n and B ∈^m × d, let = min_-n C' C' - [ A,   B ] _F, for parameters ϵ > 0 , δ > 0. The goal is to output X'∈^n× d so that there exists A' ∈^m × n  
such that

    [ A',   A'X' ] - [ A,   B ] _F ≤ (1+ϵ)  + δ.





One could solve total least squares directly, but it is much slower than solving least squares (LS). 
We will use fast randomized algorithms, the basis of which are sampling and sketching ideas <cit.>, 
to speed up solving total least squares in both theory and in practice.


The total least squares problem with regularization is also an important variant of this problem <cit.>. 
We consider the following version of the regularized total least squares problem.

Given two matrices A ∈^m × n and B ∈^m × d and λ>0, let = min_U∈^m × n, V∈^n× (n+d) UV - [ A,   B] _F^2+λU_F^2+λV_F^2, for parameters ϵ > 0 , δ > 0. The goal is to output X'∈^n× d so that there exist A' ∈^m × n  
, U'∈^m × n and V'∈^n× (n+d) satisfying [A',   A'X']-U'V'_F^2≤δ
and [ A',   A'X' ] - [ A,   B ] _F^2 ≤ (1+ϵ)  + δ.



[!t]Our Fast Total Least Squares Algorithm











§ FAST TOTAL LEAST SQUARES ALGORITHM


We present our algorithm in Algorithm <ref> and give the analysis here. 
Readers can refer to Table <ref> to check notations in Algorithm <ref>.
To clearly give the intuition,
we present a sequence of approximations, reducing the size of our problem step-by-step. We can focus on the case when d ≫Ω( n/ϵ ) and the optimal solution C to program (<ref>) has the form [ A,   AX ] ∈^m × (n+d). For the other case when d = O(n/ϵ), we do not need to use the sampling matrix D_1. In the case when the solution does not have the form [ A,   AX], we need to include Split in the algorithm, since it will perturb some columns in A with arbitrarily small noise to make sure A has rank n. By applying procedure Split, we can handle all cases.

Fix A∈^m× n and B∈^m× d.
Let = min_-n C' ∈^m × (n + d) C' - [ A,   B] _F. 
By using techniques in low-rank approximation,
we can find an approximation of a special form. More precisely, 
let S_1 ∈^s_1 × m be a CountSketch matrix with s_1 = O( n/ϵ ).
Then we claim that it is sufficient to look at solutions of the form US_1C.
[CountSketch matrix for low rank approximation problem]
With probability 0.98,

    min_-n U ∈^m × s_1 U S_1 C - C _F^2 ≤ (1+ϵ)^2 ^2.


We provide the proof in Appendix <ref>.
We shall mention that we cannot use leverage score sampling here,
because taking leverage score sampling on matrix C would take at least (C)+(n+d)^2 time,
while we are linear in d in the additive term in our running time Õ((C))+d·(n/ϵ).


Let U_1 be the optimal solution of the program min_ U ∈^m × s_1 U S_1 C - C _F^2, i.e.,

    U_1 = min_-n U ∈^m × s_1 U S_1 C - C _F^2.

If d is large compared to n,
then program (<ref>) is computationally expensive to solve.
So we can apply sketching techniques to reduce the size of the problem. 
 Let D_1^⊤∈^d_1 × (n+d) denote a leverage score sampling and rescaling matrix according to the columns of S_1 C, with d_1 = Õ(n/ϵ) nonzero entries on the diagonal of D_1.
Let U_2 ∈^m × s_1 denote the optimal solution to the problem min_-n U ∈^m × s_1 U S_1 C D_1 - C D_1 _F^2, i.e.,

    U_2 = min_-n U ∈^m × s_1 U S_1 C D_1 - C D_1 _F^2.

Then the following claim comes from the constrained low-rank approximation result (Theorem <ref>).
[Solving regression with leverage score sampling]
Let U_1 be defined in Eq. (<ref>), and let U_2 be defined in Eq. (<ref>). Then with probability 0.98,

    U_2 S_1 C - C _F^2 ≤     (1 + ϵ)^2  U_1 S_1 C - C _F^2 .


We provide the proof in Appendix <ref>. 
We now consider how to solve program (<ref>).
We observe that 

U_2 ∈colspan(C D_1).

We can thus consider the following relaxation: given C D_1, S_1 C and C, solve:

    min_-n Z ∈^d_1 × s_1 C D_1 Z S_1 C - C _F^2.

By setting CD_1Z=U,
we can check that program (<ref>) is indeed a relaxation of program (<ref>).
Let Z_1 be the optimal solution to program (<ref>). We show the following claim and delayed the proof in <ref>.
[Approximation ratio of relaxation]
With probability 0.98,

    C D_1 Z_1 S_1 C - C _F^2 ≤ (1+O(ϵ)) ^2 .


However,
program (<ref>) still has a potentially large size, i.e., we need to work with an m× d_1 matrix CD_1.
To handle this problem,
we again apply sketching techniques.
Let D_2 ∈^d_2 × m be a leverage score sampling and rescaling matrix according to the matrix C D_1 ∈^m × d_1, so that D_2 has d_2 = Õ(n/ϵ) nonzeros on the diagonal.
Now,
we arrive at the small program that we are going to directly solve: 

    min_-n Z ∈^d_1 × s_1 D_2 C D_1 Z S_1 C - D_2 C _F^2.

We shall mention that here it is beneficial to apply leverage score sampling matrix because we only need to compute leverage scores of a smaller matrix CD_1,
and computing D_2C only involves sampling a small fraction of the rows of C.
On the other hand,
if we were to use the CountSketch matrix,
then we would need to touch the whole matrix C when computing D_2C.
Overall,
using leverage score sampling at this step can reduce the constant factor of the (C) term in the running time,
and may be useful in practice.
Let -n Z_2 ∈^d_1 × s_1 be the optimal solution to this problem.

[Solving regression with a CountSketch matrix]
With probability 0.98,

    C D_1 Z_2 S_1 C - C _F^2  
    ≤  (1+ϵ)^2   C D_1 Z_1 S_1 C - C _F^2


We provide the proof in Appendix <ref>.

Our algorithm thus far is as follows: 
we compute matrices S_1, D_1, D_2 accordingly,
then solve program (<ref>) to obtain Z_2.
At this point,
we are able to obtain the low rank approximation C = CD_1 · Z_2 · S_1 C.  We show the following claim and delayed the proof in Appendix <ref>.
[Analysis of C]
With probability 0.94,

    C - C _F^2  
    ≤  (1+O(ϵ)) ^2.


Let C = [A,  B] where A∈^m × n and B∈^m × d.
However,
if our goal is to only output a matrix X so that A X =B,
then we can do this faster by not computing or storing the matrix C.
Let S_2 ∈^s_2 × m be a CountSketch matrix with s_2 = O( n/ϵ ).
 We solve a regression problem:


    min_ X ∈^n × d S_2A X - S_2B_F^2.


Notice that S_2A and S_2B are computed directly from CD_1, Z_2, S_1C and S_2.
Let X be the optimal solution to the above problem. 
[Approximation ratio guarantee]
Assume C = [A,  AX] for some X∈^n× d.
Then with probability at least 0.9,

    [A,  AX] - [A,  B] _F^2 ≤ (1+O(ϵ)) ^2.


We provide the proof in Appendix <ref>.

If the assumption C = [A,  AX] in Claim <ref> does not hold,
then we need to apply procedure Split.
Because (C)=n from our construction,
if the first n columns of C cannot span the last d columns,
then the first n columns of C are not full rank.
Hence we can keep adding a sufficiently small multiple of one of the last d columns that cannot be spanned
to the first n columns until the first n columns are full rank.
Formally,
we have

Fix s_1 =O(n/ϵ), s_2=O(n/ϵ), d_1= Õ(n/ϵ).
Given CD_1∈^m× d_1,
Z_2∈^d_1× s_1, and 
S_1C∈^s_1× (n+d)
so that C:=CD_1· Z_2· S_1C has rank n,
procedure SPLIT (Algorithm <ref>) returns A∈^s_2× n and B∈^s_2× d in time O((C)+d·(n/ϵ)) so that there exists X∈^n× d satisfying
A·X=B.
Moreover,
letting A be the matrix computed in lines (<ref>) to (<ref>),
then with probability 0.99,

    [A,  AX]-C_F≤C-C_F+δ.


We provide the proof in Appendix <ref>. 
Now that we have A and B, and we can compute X by solving the regression problem min_X∈^n× dA X - B_F^2.

We next summarize the running time. Ommitted proofs are in Appendix <ref>. 

Procedure FastTotalLeastSquares in Algorithm <ref>  runs in time Õ((A) + (B) + d ·(n/ϵ)).

To summarize, Theorem <ref> shows the performance of our algorithm. Ommitted proofs are in Appendix <ref>.

Given two matrices A ∈^m × n and B ∈^m × d, letting 

    = min_-n C' ∈^m × (n + d) C' - [ A,   B] _F,

we have that for any ϵ∈ (0,1), there is an algorithm (procedure FastTotalLeastSquares in Algorithm <ref>) that runs in Õ((A) + (B)) + d ·(n/ϵ) time and outputs a matrix X ∈^n × d such that there is a matrix A∈^m × n satisfying that

    [A,  AX] - [A,   B] _F ≤ (1+ϵ)  + δ

holds with probability at least 9/10, where δ > 0 is arbitrarily small.



The success probability 9/10 in Theorem <ref> can be boosted to 1-δ for any δ>0 in a standard way. Namely, we run our FTLS algorithm O(log(1/δ)) times where in each run we use independent randomness, and choose the solution found with the smallest cost. Note that for any fixed output X, the cost [A̅, A̅ X ] - [A, B] _F can be efficiently approximated. To see this, let S be a CountSketch matrix with O(ϵ^-2) rows. Then S[A̅, A̅ X ] - S[A, B]_F = (1 ±ϵ) [A̅, A̅ X ] - [A, B] _F with probability 9/10 (see, for example Lemma 40 of <cit.> ). We can compute S[A̅, A̅ X ] - S[A, B]_F in time O(d ·poly(n/ϵ)), and applying S can be done in nnz(A) + nnz(B) time. We can then amplify the success probability by taking O(log (1/δ)) independent estimates and taking the median of the estimates. This is a (1±ϵ)-approximation with probability at least 1- O(δ/log(1/δ)). We run our FTLS algorithm O(log(1/δ)) times, obtaining outputs X^1, …, X^O(log(1/δ)) and for each X^i, apply the method above to estimate its cost. Since for each X^i our estimate to the cost is within 1 ±ϵ with probability at least 1-O(δ/(log(1/δ)), by a union bound the estimates for all X^i are within 1 ±ϵ with probability at least 1-δ/2. Since also the solution with minimal cost is a 1 ±ϵ approximation with probability at least 1-δ/2, by a union bound we can achieve 1-δ probability with running time Õ(log^2 (1/δ)) · (nnz(A) + nnz(B) + d·poly(n/ϵ))).



We further generalize our algorithm to handle regularization. Ommitted proofs can be found in Appendix <ref>.

Given two matrices A ∈^m × n and B ∈^m × d and λ>0, letting 

    = min_U∈^m × n, V∈^n× (n+d) UV - [ A,   B] _F^2+λU_F^2+λV_F^2,
 
we have that for any ϵ∈ (0,1), there is an algorithm (procedure FastRegularizedTotalLeastSquares in Algorithm <ref>) 
 that runs in 

    Õ((A) + (B) + d ·(n/ϵ))

time and outputs a matrix X ∈^n × d such that there is a matrix A∈^m × n, U∈^m × n and V∈^n× (n+d) satisfying [A,  AX]-UV_F^2≤δ and with probability 9/10, 

    [A,  AX] - [A,   B] _F^2 + λU_F^2 + λV_F^2 ≤ (1+ϵ)  + δ.












§ EXPERIMENTS



We conduct several experiments to verify the running time and optimality of our fast total least squares algorithm <ref>.
Let us first recall the multiple-response regression problem.
Let A∈ℝ^m× n and B∈ℝ^m× d.
In this problem, we want to find X∈ℝ^n× d so that AX∼ B.
The least squares method (LS) solves the following optimization program:

    c_:=    min_ X∈ℝ^n× d, Δ B∈ℝ^m× dΔ B _F^2, 
       subject to   AX=B+Δ B.

On the other hand, the total least squares method (TLS) solves the following optimization program:

    c_:=min_-n C' ∈^m × (n+d) C' - [A   B] _F.

The fast total least squares method (FTLS) returns X∈^n× d,
which provides an approximation C'=[A AX] to the TLS solution,
and the cost is computed as c_ =   C' - C _F^2.

Our numerical tests are carried out on an Intel Xeon E7-8850 v2 server with 2.30GHz and 4GB RAM under Matlab R2017b.
[The code can be found at <https://github.com/yangxinuw/total_least_squares_code>.]



 §.§ A Toy Example


We first run our FTLS algorithm on the following toy example, for which we have the analytical solution exactly.
Let A∈ℝ^3× 2 be A_11=A_22=1 and 0 everywhere else.
Let B∈ℝ^3× 1 be B_3=3 and 0 everywhere else.
We also consider the generalization of this example with larger dimension in Appendix <ref>.
The cost of LS is 9,
since AX can only have non-zero entries on the first 2 coordinates,
so the 3rd coordinate of AX-B must have absolute value 3. Hence the cost is at least 9.
Moreover,
a cost 9 can be achieved by setting X=0 and Δ B=-B.
However,
for the TLS algorithm,
the cost is only 1.
Consider Δ A∈ℝ^3× 2 where A_11=-1 and 0 everywhere else.
Then C':=[(A+Δ A),   B] has rank 2,
and C'-C_F=1.


We first run experiments on this small matrix.
Since we know the solution of LS and TLS exactly in this case,
it is convenient for us to compare their results with that of the FTLS algorithm.
When we run the FTLS algorithm,
we sample 2 rows in each of the sketching algorithms.

The experimental solution of LS is C_=(0,1,3) which matches the theoretical solution. The cost is 9. The experimental solution of TLS is C_= (1,1,0) which also matches the theoretical result.
The cost is 1. 

FTLS is a randomized algorithm,
so the output varies.
We post several outputs:

    C_=
    [  .06 -.01  .25; -.01  .99  .00;  .76  .01 2.79;      ],
    [  .14 -.26 -.22; -.26  .91 -.06; -.67 -.20 2.82;      ]
 
These solutions have cost of 1.55 and  1.47.





    Ĉ_=
    [  0.2370  0.3470 -0.2458;  0.3470  0.8422  0.1118; -0.7373  0.3353  2.7625;         ]

This solution has a cost of 1.6333.























We run the FTLS multiple times to analyze the distribution of costs.
Experimental result,
which can be found in Appendex <ref>,
shows that FTLS is a stable algorithm,
and consistently performs better than LS. 

We also consider a second small toy example.
Let A still be a 10× 5 matrix and B be a 10× 1 vector.
Each entry A(i,j) is chosen i.i.d. from the normal distribution N(0,1),
and each entry B(i) is chosen from N(0,3). 
Because entries from A and B have different variance,
we expect the results of LS and TLS to be quite different.
When we run the FTLS algorithm,
we sample 6 rows.

We run FTLS 1000 times, and compute the distribution of costs.
The results of this experiment,
which is in Appendex <ref>,
again demonstrates the stability of the algorithm.



 §.§ Large Scale Problems







We have already seen that FTLS works pretty well on small matrices.
We next show that the fast total least squares method also provides a good estimate 
for large scale regression problems.
The setting for matrices is as follows:
for k=5,10,⋯,100,
we set A to be a 20k× 2k matrix where A(i,i)=1 for i=1,⋯,2k and 0 everywhere else, 
and we set B to be a 20k× 1 vector where B(2k+1)=3 and 0 elsewhere.
As in the small case,
the cost of TLS is 1,
and the cost of LS is 9.

Recall that in the FTLS algorithm,
we use Count-Sketch/leverage scores sampling/Gaussian sketches to speed up the algorithm.
In the experiments,
we take sample density ρ=0.1,0.3,0.6,0.9 respectively to check our performance.
The left 2 pictures in Figure <ref> show 
the running time together with the ratio TLS/FTLS for different sample densities.

We can see that the running time of FTLS is significantly smaller than that of TLS.
This is because the running time of TLS depends heavily on m,
the size of matrix A.
When we apply sketching techniques,
we significantly improve our running time.
The fewer rows we sample,
the faster the algorithm runs.
We can see that FTLS has pretty good performance;
even with 10% sample density,
FTLS still performs better than LS.
Moreover,
the more we sample,
the better accuracy we achieve.





















The above matrix is extremely sparse.
We also consider another class of matrices.
For k=5,10,⋯,100,
we set A to be a 20k× 2k matrix where A(i,j)∼ N(0,1);
we set B to be a 20k× 1 vector where B(i)∼ N(0,3). 
As in previous experiments,
we take sample densities of ρ=0.1,0.3,0.6,0.9, respectively, to check our performance.
The results of this experiment are shown in the right 2 pictures in Figure <ref>.

We see that compared to TLS,
our FTLS sketching-based algorithm significantly reduces the running time.
FTLS is still slower than LS,
though,
because in the FTLS algorithm we still need to solve a LS problem of the same size.
However, as discussed, LS is inadequate in a number of applications 
as it does not allow for changing the matrix A. 


The accuracy of our FTLS algorithms is also shown.














We also conducted experiments on real datasets from the UCI Machine Learning Repository <cit.>.
We choose datasets with regression task.
Each dataset consists of input data and output data.
To turn it into a total least squares problem,
we simply write down the input data as a matrix A and the output data as a matrix B,
then run the corresponding algorithm on (A,B). 
We have four real datasets :
 Airfoil Self-Noise <cit.> in Table <ref>(a),
 Wine Quality Red wine <cit.> in Table <ref>(b),
 Wine Quality White wine <cit.> in Table <ref>(c),
 Insurance Company Benchmark (COIL 2000) Data Set <cit.>
 From the results,,
 we see that FTLS also performs well on real data:
 when FTLS samples 10% of the rows,
 the result is within 5% of the optimal result of TLS,
 while the running time is 20-40 times faster.
 In this sense,
 FTLS achieves the advantages of both TLS and LS: FTLS has almost the same accuracy as TLS,
 while FTLS is significantly faster. 





















alpha







§ APPENDIX




§ NOTATION

In addition to O(·) notation, for two functions f,g, we use the shorthand f≲ g (resp. ≳) to indicate that f≤ C g (resp. ≥) for an absolute constant C. We use f g to mean cf≤ g≤ Cf for constants c,C.



§ OBLIVIOUS AND NON-OBLIVIOUS SKETCHING MATRIX


In this section we introduce techniques in sketching.
In order to optimize performance,
we introduce multiple types of sketching matrices,
which are used in Section <ref>.
In Section <ref>, we provide the definition of CountSketch and Gaussian Transforms. In Section <ref>, we introduce leverage scores and sampling based on leverage scores.



 §.§ CountSketch and Gaussian Transforms


CountSketch matrix comes from the data stream literature <cit.>.

A CountSketch transform is defined to be Π=  Φ D∈ℝ^m× n. Here,  D is an n× n random diagonal matrix with each diagonal entry independently chosen to be +1 or -1 with equal probability, and Φ∈{0,1}^m× n is an m× n binary matrix with Φ_h(i),i=1 and all remaining entries 0, where h:[n]→ [m] is a random map such that for each i∈ [n], h(i) = j with probability 1/m for each j ∈ [m]. For any matrix A∈ℝ^n× d, Π A can be computed in O((A)) time. 


To obtain the optimal number of rows,
we need to apply Gaussian matrix,
which is another well-known oblivious sketching matrix.

Let S=1/√(m)· G ∈ℝ^m× n where  each entry of G∈ℝ^m× n is chosen independently from the standard Gaussian distribution. For any matrix A∈ℝ^n× d, SA can be computed in O(m ·(A)) time. 


We can combine CountSketch and Gaussian
transforms to achieve the following:

Let S' = S Π, where Π∈ℝ^t× n is the CountSketch transform (defined in Definition <ref>) and S∈ℝ^m × t is the Gaussian transform (defined in Definition <ref>). For any matrix A∈ℝ^n× d, S'A can be computed in O((A) + dtm^ω-2) time, where ω is the matrix multiplication exponent.




 §.§ Leverage Scores


We do want to note that there are other ways of constructing sketching matrix though, such as through sampling the rows of A via
a certain distribution and reweighting them.
This is called leverage score sampling <cit.>. 
We first give the concrete definition of leverage scores.

Let U∈ℝ^n× k have orthonormal columns with n ≥ k. We will use the notation p_i = u_i^2 / k, where u_i^2 =  e_i^⊤ U _2^2 is referred to as the i-th leverage score of U. 



Next we explain the leverage score sampling.
 Given A∈ℝ^n× d with rank k, let U∈ℝ^n× k be an orthonormal basis of the column span of A, and for each i let k· p_i be the squared row norm of the i-th row of U. Let p_i denote the i-th leverage score of U. Let β>0 be a constant and q=(q_1, ⋯,q_n) denote a distribution such that, for each i∈ [n], q_i ≥β p_i. Let s be a parameter. Construct an n× s sampling matrix B and an s× s rescaling matrix D as follows. Initially, B = 0^n× s and D = 0^s× s. For the same column index j of B and of D, independently, and with replacement, pick a row index i∈ [n] with probability q_i, and set B_i,j=1 and D_j,j=1/√(q_i s). We denote this procedure Leverage score sampling according to the matrix A.


Leverage score sampling is efficient in the sense that leverage score can be efficiently approximated.

For any ϵ>0,
with probability at least 2/3,
we can compute 1±ϵ approximation of all leverage scores of matrix A∈^n× d in time O((A)+r^ωϵ^-2ω) where r is the rank of A and ω≈ 2.373 is the exponent of matrix multiplication <cit.>.



In Section  <ref> we show how to apply matrix sketching to solve regression problems faster.
In Section <ref>,
we give a structural result on rank-constrained approximation problems.


§ MULTIPLE REGRESSION

Linear regression is a fundamental problem in Machine Learning. There are a lot of attempts trying to speed up the running time of different kind of linear regression problems via sketching matrices <cit.>. A natural generalization of linear regression is multiple regression.

We first show how to use CountSketch to reduce to a multiple regression problem:

Given A ∈^n × d and B ∈^n × m, let S ∈^s × n denote a sampling and rescaling matrix according to A. Let X^* denote min_X  A X - B _F^2 and X' denote min_X  S A X - S B _F^2.  If S has s = O(d / ϵ) rows, then we have that

    A X' - B _F^2 ≤ (1+ϵ)  A X^* - B _F^2

holds with probability at least 0.999.


The following theorem says leverage score sampling solves multiple response 
regression:

Given A ∈^n × d and B ∈^n × m, let D ∈^n × n denote a sampling and rescaling matrix according to A. Let X^* denote min_X  A X - B _F^2 and X' denote min_X  D A X - S B _F^2.  If D has O(dlog d + d/ϵ) non-zeros in expectation, that is, this is the
expected number of sampled rows, then we have that

    A X' - B _F^2 ≤ (1+ϵ)  A X^* - B _F^2

holds with probability at least 0.999.





§ GENERALIZED RANK-CONSTRAINED MATRIX APPROXIMATION


We state a tool which has been used in several recent works <cit.>.

Given matrices A ∈^n × d, B ∈^n × p, and C ∈^q × d, let the singular value decomposition (SVD) of B be B = U_B Σ_B V_B^⊤ and the SVD of C be C = U_C Σ_C V_C^⊤. Then

    B^† ( U_B U_B^⊤ A V_C V_C^⊤ )_k C^† = -k X ∈^p × qmin A - B X C _F

where (U_B U_B^⊤ A V_C V_C^⊤)_k ∈^n × d is of rank at most k and denotes the best -k approximation to U_B U_B^⊤ A V_C V_C^⊤∈^n × d in Frobenius norm.

Moreover,
( U_B U_B^⊤ A V_C V_C^⊤ )_k can be computed by first computing the SVD decomposition of U_B U_B^⊤ A V_C V_C^⊤ in time O(nd^2),
then only keeping the largest k coordinates.
Hence B^† ( U_B U_B^⊤ A V_C V_C^⊤ )_k C^† can be computed in O(nd^2+np^2+qd^2) time.





§ CLOSED FORM FOR THE TOTAL LEAST SQUARES PROBLEM


Markovsky and Huffel <cit.> propose the following alternative formulation of total least squares problem.

    min_-n C' ∈^m × (n+d) C' - C _F

When program (<ref>) has a solution (X,Δ A,Δ B),
we can see that (<ref>) and (<ref>) are in general equivalent by setting C'=[A + Δ A,    B + Δ B].
However,
there are cases when program (<ref>) fails to have a solution,
while (<ref>) always has a solution.


As discussed, a solution to the total least squares problem
can sometimes be written in closed form.  
Letting C=[A,   B],
denote the singular value decomposition (SVD) of C by UΣ V^⊤, where Σ=𝖽𝗂𝖺𝗀(σ_1,⋯,σ_n+d)∈^m× (n+d) with σ_1≥σ_2≥⋯≥σ_n+d. 
Also we represent (n+d)× (n+d) matrix V as [ V_11 V_12; V_21 V_22;      ]
where V_11∈^n× n and V_22∈^d× d.

Clearly Ĉ=U𝖽𝗂𝖺𝗀(σ_1,⋯,σ_n,0,⋯,0)V^⊤ is a minimizer of program (<ref>). 
But whether a solution to program (<ref>) exists depends on the singularity of V_22.
In the rest of this section we introduce different cases of the solution to program (<ref>), and discuss how our algorithm deals with each case.





 §.§ Unique Solution

We first consider the case when 
the Total Least Squares problem has a unique solution.

If σ_n>σ_n+1, and
V_22 is non-singular,
then the minimizer Ĉ is given by U𝖽𝗂𝖺𝗀(σ_1,⋯,σ_n,0,⋯,0)V^⊤,
and the optimal solution X̂ is given by -V_12V_22^-1.


Our algorithm will first find a rank n matrix C'=[A',   B'] so that C'-C_F is small,
then solve a regression problem to find X' so that A'X'= B'.
In this sense,
this is the most favorable case to work with,
because a unique optimal solution Ĉ exists,
so if C' approximates Ĉ well,
then the regression problem A'X'=B' is solvable.



 §.§ Solution exists, but is not unique


If σ_n=σ_n+1,
then it is still possible that the Total Least Squares problem has a unique solution,
although this time, the solution X̂ is not unique.
Theorem <ref> is a generalization of Theorem <ref>.

Let p≤ n be a number so that σ_p>σ_p+1=⋯ =σ_n+1.
Let V_p be the submatrix that contains the last d rows and the last n-p+d columns of V. 
If 
V_p is non-singular,
then multiple minimizers Ĉ=[Â,  B̂] exist,
and there exists X̂∈^n× d so that ÂX̂=B̂.


We can also handle this case. 
As long as the Total Least Squares problem has a solution X̂,
we are able to approximate it by first finding C'=[A',   B'] and then solving a regression problem. 



 §.§ Solution does not exist

Notice that the cost C-C_F^2, where C is the optimal solution to program (<ref>),
always lower bounds the cost of program (<ref>).
But there are cases where this cost is not approchable in program (<ref>).

If V_22 is singular,
letting C denote [A,  B],
then AX=B has no solution.

Theorem <ref> shows that even if we can compute Ĉ precisely,
we cannot output X,
because the first n columns of Ĉ cannot span the rest d columns.
In order to generate a meaningful result,
our algorithm will perturb C' by an arbitrarily small amount 
so that A'X'=B' has a solution.
This will introduce an arbitrarily small additive error  
in addition to our relative error guarantee. 























































§ OMITTED PROOFS IN SECTION <REF>





 §.§ Proof of Claim <ref>



Let C^* be the optimal solution of min_-n C' ∈^m × (n + d) C' - [ A,   B] _F.
Since (C^*)=n≪ m,
there exist U^*∈^m× s_1 and V^*∈^s_1× (n+d) so that C^*=U^*V^*,
and (U^*)=(V^*)=n.

Therefore

    min_V∈^s_1× (n+d) U^* V - C _F^2 =^2.

Now consider the problem formed by multiplying by S_1 on the left,

    min_V∈^s_1× (n+d) S_1 U^* V - S_1 C _F^2 .

Letting V' be the minimizer to the above problem, we have

    V' = (S_1 U^*)^† S_1 C.


Thus, we have

    min_-n U∈^m × s_1 U S_1 C - C _F^2 
    ≤     U^* (S_1 U^*)^† S_1 C - C _F^2 
    
    =      U^* V' - C _F^2 
    ≤     (1+ϵ)  S_1 U^* V' - S_1 C _F^2 
    ≤     (1+ϵ)  S_1 U^* V^* - S_1 C _F^2 
    ≤     (1+ϵ)^2  U^* V^* - C _F^2 
    
    =      (1+ϵ)^2 ^2

where the first step uses the fact that U^* (S_1 U^*)^† S_1∈^m × s_1 with rank n,
the second step is the definition of V',
the third step follows from the definition of the Count-Sketch matrix S_1 and Theorem <ref>,
the fourth step uses the optimality of V',
and the fifth step again uses Theorem <ref>.





 §.§ Proof of Claim <ref>




We have

    U_2 S_1 C - C _F^2 ≤     (1+ϵ)  U_2 S_1 C D_1- C D_1_F^2 
    ≤     (1+ϵ)  U_1 S_1 C D_1- C D_1_F^2 
    ≤     (1 + ϵ)^2  U_1 S_1 C - C _F^2,

where the first step uses the property of a leverage score sampling matrix D_1, the second step follows from the definition of U_2 (i.e., U_2 is the minimizer), and the last step follows from the property of the leverage score sampling matrix D_1 again. 





 §.§ Proof of Claim <ref>




From Claim <ref> we have 
that U_2 ∈colspan(C D_1). 
Hence we can choose Z so that CD_1Z=U_2.
Then by Claim <ref> and Claim <ref>,
we have

    C D_1 Z S_1 C - C _F^2=  U_2 S_1 C - C _F^2 ≤ (1+ϵ)^4 ^2.

Since Z_1 is the optimal solution, the objective value can only be smaller.




 §.§ Proof of Claim <ref>




Recall that Z_1=min_-n Z ∈^d_1 × s_1 C D_1 Z S_1 C - C _F^2.
Then we have

    C D_1 Z_2 S_1 C - C _F^2  
    ≤    (1+ϵ)  D_2C D_1 Z_2 S_1 C - D_2C _F^2
    ≤    (1+ϵ)  D_2C D_1 Z_1 S_1 C - D_2C _F^2
    ≤     (1+ϵ)^2  C D_1 Z_1 S_1 C - C _F^2,

where the first step uses the property of the leverage score sampling matrix D_2,
the second step follows from the definition of Z_2 (i.e., Z_2 is a minimizer), 
and the last step follows from the property of the leverage score sampling matrix D_2.






 §.§ Proof of Claim <ref>




    C - C _F^2  
     =     CD_1 · Z_2 · S_1 C - C _F^2
    ≤    (1+ϵ)^2  C D_1 Z_1 S_1 C - C _F^2 
    ≤     (1+O(ϵ)) ^2

where the first step is the definition of C,
the second step is Claim <ref>,
and the last step is Claim <ref>.




 §.§ Proof of Claim <ref>



By the condition that C=[A AX],
B=AX,
hence
X is the optimal solution to the program min_ X ∈^n × dA X - B_F^2.
Hence by Theorem <ref>,
with probability at least 0.99,

    AX - B_F^2≤ (1+ϵ)AX - B_F^2=0

Therefore

    [A,  AX] - [A,  B] _F^2 =  [A,  B ] - C _F^2= C  - C _F^2.

Then it follows from Claim <ref>.




 §.§ Proof of Lemma <ref>




Proof of running time.
Let us first check the running time.
We can compute C=S_2·C by first computing S_2· CD_1,
then computing (S_2CD_1)· Z_2,
then finally computing S_2CD_1Z_2S_1C.
Notice that D_1 is a leverage score sampling matrix,
so (CD_1)≤(C).
So by Definition <ref>,
we can compute S_2· CD_1 in time O((C)).
All the other matrices have smaller size,
so we can do matrix multiplication in time O(d·(n/ϵ)).
Once we have C,
the independence between columns in A can be checked in time O(s_2· n).
The FOR loop will be executed at most n times,
and inside each loop,
line (<ref>) will take at most d linear independence checks.
So the running time of the FOR loop is at most O(s_2· n)· n· d=O(d·(n/ϵ)).
Therefore the running time is as desired. 

Proof of Correctness. 
We next argue the correctness of procedure Split.
Since (C)=n,
with high probability (C)=(S_2·C)= n.
Notice that B is never changed in this subroutine.
In order to show there exists an X so that AX=B,
it is sufficient to show that at the end of procedure Split,
(A)=(C),
because this means that the columns of A span each of the columns of C,
including B.
Indeed,
whenever (A_*,[i])<i,
line 25 will be executed.
Then by doing line 26,
the rank of A will increase by 1,
since by the choice of j,
A_*,i+δ·B_*,j is independent form A_*,[i-1].
Because (C)=n,
at the end of the FOR loop we will have (A)=n.

Finally let us compute the cost.
In line (<ref>) we use δ/(m),
and thus 

    [A,  B]- C_F^2≤δ^2/(m)·B_F^2≤δ^2.

We know that
X is the optimal solution to the program min_ X ∈^n × d S_2A X - S_2B_F^2.
Hence by Theorem <ref>,
with probability 0.99,

    AX - B_F^2≤ (1+ϵ)min_ X ∈^n × d S_2A X - S_2B_F^2=0.

which implies AX= B.
Hence we have

    [A,  AX] - C _F  
    ≤     [A,  AX] - C_F+C - C _F
    
    =      [A,  B] - C_F+C - C _F
    ≤    δ+C - C _F

where the first step follows by triangle inequality, and the last step follows by (<ref>).
 




 §.§ Proof of Lemma <ref>




We bound the time of each step:

1. Construct the s_1× m Count-Sketch matrix S_1 and compute S_1C with s_1=O( n/ϵ ). This step takes time (C)  + d ·(n/ϵ).

2. Construct the (n+d) × d_1 leverage sampling and rescaling matrix D_1  with d_1 = Õ(n/ϵ) nonzero diagonal entries and compute CD_1. This step takes time Õ((C)  + d ·(n/ϵ)).

3. Construct the d_2 × m leverage sampling and rescaling matrix D_2  with d_2 = Õ(n/ϵ) nonzero diagonal entries. This step takes time Õ((C)  + d ·(n/ϵ)) according to Theorem <ref>.

4. Compute Z_2 ∈^d_1 × s_1 by solving the rank-constrained system: 

    min_-n Z ∈^d_1 × s_1 D_2 C D_1 Z S_1 C - D_2 C _F^2.

Note that D_2 C D_1 has size Õ(n/ϵ)×Õ(n/ϵ), S_1C has size O( n/ϵ )× (n+d), and D_2C has size Õ(n/ϵ)× (n+d), so according to Theorem <ref>, we have an explicit closed form for Z_2, and the time taken is d·(n/ϵ).

5. Run procedure Split to get A∈^s_2× n and B∈^s_2× d with s_2=O( n/ϵ ).
By Lemma <ref>,
this step takes time O((C)+d·(n/ϵ)).

6. Compute X by solving the regression problem min_X∈^n× dA X - B_F^2 in time O(d ·(n/ϵ)). This is because X= ( A )^†B, and A has size O( n/ϵ )× n, so we can compute ( A )^† in time O(( n/ϵ )^ω)=(n/ϵ), and then compute X in time O( (n/ϵ)^2· d ) since B is an O( n/ϵ )× d matrix.


	



Notice that (C)=(A)+(B), so we have the desired running time. 




 §.§ Procedure Evaluate


In this subsection we explain what procedure Evaluate does.
Ideally,
we would like to apply procedure Split on the matrix C directly so that the linear system 
AX=B has a solution.
However,
C has m rows, which is computationally expensive to work with.
So in the main algorithm we actually apply procedure Split on the sketched matrix S_2C. 
When we need to compute the cost,
we shall redo the operations in procedure Split on C to split C correctly.
This is precisely what we are doing in lines (<ref>) to (<ref>).



 §.§ Putting it all together


The running time follows from Lemma <ref>. For the approximation ratio, 
let A, A be defined as in Lemma <ref>. 
From Lemma <ref>, there exists X∈^n × d satisfying AX = B. Since X is obtained from solving the regression problem A X - B_F^2, we also have A X = B. Hence with probability 0.9, 

[ A, A X ] - C _F 
≤δ + C - C _F 
≤δ + ( 1+ O ( ϵ ) ) , 

where the first step uses Lemma <ref> and the second step uses Claim <ref>.
Rescaling ϵ gives the desired statement.






§ EXTENSION TO REGULARIZED TOTAL LEAST SQUARES PROBLEM



In this section we provide our algorithm for the regularized total least squares problem and prove its correctness.
Recall our regularized total least squares problem is defined as follows.

    :=min_A∈^m× n, X∈^n× d,U∈^m × n, V∈^n× (n+d)    UV - [ A,   B] _F^2+λU_F^2+λV_F^2
     subject to     [A, AX]=UV

[!t]Our Fast Total Least Squares Algorithm with Regularization




For λ>0 and rank k matrix A,
the statistical dimension of the ridge regression problem with regularizing weight λ is defined as

    _λ(A):=∑_i∈ [k]1/1+λ/σ_i^2

where σ_i is the i-th singular value of A for i∈[k].

Notice that _λ(A) is decreasing in λ,
so we always have _λ(A)≤_0(A)=(A).


Given positive integers n_1, n_2, r, s, k and parameter λ≥ 0. 
For C ∈^n_1 × r, D ∈^s × n_2, B ∈^n_1 × n_2, the problem of finding

    min_Z_R ∈^r × k , Z_S ∈^k × s C Z_R Z_S D - B _F^2 + λ C Z_R _F^2 + λ Z_S D _F^2,

and the minimizing of C Z_R ∈^n_1 × k and Z_S D ∈^k × n_2, can be solved in

    O(n_1 r ·(C) + n_2 s ·(D) + (D) · n_1 (n_2 + r_C) )

time.












Fix m≥ n.
For A∈^m× n, B∈^n× d and λ>0,
consider the rigid regression problem

    min_X∈^n× dAX-B_F^2+λX_F^2.

Let S ∈^s × m be a CountSketch matrix with s = O( _λ(A)/ϵ )=O(n/ϵ ),
then with probability 0.99,

    min_X∈^n× dSAX-SB_F^2+λX_F^2≤ (1+ϵ)min_X∈^n× dAX-B_F^2+λX_F^2

Moreover,
SA, SB can be computed in time 

    O((A)+(B))+O((n+d)(_λ(A)/ϵ+_λ(A)^2)).




We claim that it is sufficient to look at solutions of the form CS_2^⊤ Z_2Z_1S_1C.
[CountSketch matrix for low rank approximation problem]
Given matrix C ∈^m × (n+d). Let  be defined as in (<ref>). For any ϵ>0,
let S_1∈^s_1× m, S_2∈^s_2× m be the sketching matrices defined in Algorithm <ref>,
then with probability 0.98,

    min_ Z_1∈^n× s_1,Z_2∈^s_2× n CS_2^⊤ Z_2Z_1S_1C - C _F^2+λCS_2^⊤ Z_2_F^2+λZ_1S_1C_F^2 ≤ (1+ϵ)^2 .



Let U^*∈^m× n and V^*∈^n× (n+d) be the optimal solution to the program (<ref>).
Consider the following optimization problem:

    min_ V∈^n× (n+d) U^*V - C _F^2+λV_F^2


Clearly V^* ∈^n × (n+d) is the optimal solution to program (<ref>),
since for any solution V ∈^n × (n+d) to program (<ref>) with cost c,
(U^*,V) is a solution to program (<ref>) with cost c+λU^*_F^2.

Program (<ref>) is a ridge regression problem.








Hence we can take a CountSketch matrix S∈^s_1× m with s_1=O(n/ϵ) to obtain

    min_ V∈^n× (n+d) S_1 U^*V - S_1 C _F^2+λV_F^2

Let V_1 ∈^n× (n+d) be the minimizer of the above program,
then we know 

    V_1=[   S_1U^*; √(λ) I_n ]^†[ S_1C;    0 ],

which means V_1 ∈^n × (n+d) lies in the row span of S_1C ∈^ s_1 × (n+d).
Moreover,
by Theorem <ref>, with probability at least 0.99 we have

    U^*V_1 - C _F^2+λV_1_F^2≤ (1+ϵ) U^*V^* - C _F^2+λV^*_F^2











Now consider the problem

    min_U∈^m × n UV_1 - C _F^2+λU_F^2

Let U_0 ∈^m × n be the minimizer of program (<ref>).
Similarly,









 we can take a CountSketch matrix S_2∈^s_2× (n+d) with s_2=O(n/ϵ) to obtain

    min_ U∈^m× n UV_1S_2^⊤ - C S_2^⊤_F^2+λU_F^2

Let U_1 ∈^m × n be the minimizer of program (<ref>),
then we know 

    U_1^⊤=[ S_2V_1^⊤; √(λ) I_n ]^†[ S_2C^⊤;      0 ],


which means U_1 ∈^m × n lies in the column span of CS_2^⊤∈^m × s_2.
Moreover, with probability at least 0.99 we have

    U_1V_1- C  _F^2+λU_1_F^2
       ≤ (1+λ)·( U_0V_1- C  _F^2+λU_0_F^2)
       ≤ (1+λ)·( U^*V_1- C  _F^2+λU^*_F^2)

where the first step we use Theorem <ref> and the second step follows that U_0 is the minimizer.










Now let us compute the cost.

    U_1V_1 - C_F^2+λU_1_F^2+λV_1_F^2
    
    =    λV_1_F^2+( U_1V_1 - C_F^2+λU_1_F^2)
    ≤    λV_1_F^2+(1+ϵ)( U^*V_1 - C_F^2+λU^*_F^2)
    ≤    (1+ϵ)·(λU^*_F^2+( U^*V_1 - C_F^2+λV_1_F^2))
    ≤    (1+ϵ)·(λU^*_F^2+(1+ϵ)^2·( U^*V^* - C_F^2+λV^*_F^2))
    ≤    (1+ϵ)^2·( U^*V^* - C_F^2+λU^*_F^2+λV^*_F^2)
    
    =    (1+ϵ)^2

where the second step follows from (<ref>),
the fourth step follows from (<ref>),
and the last step follows from the definition of U^* ∈^m × n,V^* ∈^n × (n+d).

Finally,
since V_1 ∈^n × (n+d) lies in the row span of S_1C ∈^s_1 × (n+d) and U_1 ∈^m × n lies in the column span of CS_2^⊤∈^m × s_2,
there exists Z_1^*∈^n× s_1 and Z_2^*∈^s_2× n so that V_1=Z_1^*S_1C ∈^n × (n+d) and U_1=CS_2^⊤ Z_2^* ∈^m × n.
Then the claim stated just follows from (Z_1^*,Z_2^*) are also feasible.


Now we just need to solve the optimization problem

    min_ Z_1∈^n× s_1,Z_2∈^s_2× n CS_2^⊤ Z_2Z_1S_1C - C _F^2+λCS_2^⊤ Z_2_F^2+λZ_1S_1C_F^2

The size of this program is quite huge,
i.e., we need to work with an m× d_2 matrix CS_2^⊤.
To handle this problem,
we again apply sketching techniques.
Let D_1 ∈^d_1 × m be a leverage score sampling and rescaling matrix according to the matrix C S_2 ∈^m × s_2, so that D_1 has d_1 = Õ(n/ϵ) nonzeros on the diagonal.
Now,
we arrive at the small program that we are going to directly solve: 

    min_ Z_1∈^n× s_1,Z_2∈^s_2× n D_1CS_2^⊤ Z_2Z_1S_1C - D_1C _F^2+λD_1CS_2^⊤ Z_2_F^2+λZ_1S_1C_F^2


We have the following approximation guarantee.


Let (Z_1^*,Z_2^*) be the optimal solution to program (<ref>).
Let (Z_1,Z_2) be the optimal solution to program (<ref>).
With probability 0.96,

    CS_2^⊤Z_2Z_1S_1C - C _F^2+λCS_2^⊤Z_2_F^2+λZ_1S_1C_F^2 
    ≤   (1+ϵ)^2( CS_2^⊤ Z_2^*Z_1^*S_1C - C _F^2+λCS_2^⊤ Z_2^*_F^2+λZ_1^*S_1C_F^2)



This is because

    CS_2^⊤Z_2Z_1S_1C - C _F^2+λCS_2^⊤Z_2_F^2+λZ_1S_1C_F^2 
    ≤   (1+ϵ)( D_1CS_2^⊤Z_2Z_1S_1C - D_1C _F^2+λD_1CS_2^⊤Z_2_F^2)+λZ_1S_1C_F^2 
    ≤   (1+ϵ)( D_1CS_2^⊤Z_2Z_1S_1C - D_1C _F^2+λD_1CS_2^⊤Z_2_F^2+λZ_1S_1C_F^2) 
    ≤   (1+ϵ)( D_1CS_2^⊤ Z_2^*Z_1^*S_1C - D_1C _F^2+λD_1CS_2^⊤ Z_2^*_F^2+λZ_1^*S_1C_F^2) 
    ≤   (1+ϵ)^2( CS_2^⊤ Z_2^*Z_1^*S_1C - C _F^2+λCS_2^⊤ Z_2^*_F^2+λZ_1^*S_1C_F^2)

where the first step  uses property of the leverage score sampling matrix D_1,
the third step follows from (Z_1,Z_2) are minimizers of program (<ref>),
and the fourth step again uses property of the leverage score sampling matrix D_1.


Let U=CS_2^⊤Z_2,V=Z_1S_1C and C=UV.
Combining Claim <ref> and Claim <ref> together,
we get with probability at least 0.91,

    UV - [ A,   B] _F^2+λU_F^2+λV_F^2≤ (1+ϵ)^4

If the first n columns of C can span the whole matrix C,
then we are in good shape. 
In this case we have:
[Perfect first n columns]
Let S_3∈^s_3× m be the CountSketch matrix defined in Algorithm <ref>.
Write C as [A, B] where A∈^m× n and B∈^m× d.
If there exists X∈^n× d so that B=AX,
let X̅∈^n× d be the minimizer of min_X∈^n× dS_3AX-S_3B_F^2,
then with probability 0.9,

    [A, AX̅] - [ A,   B] _F^2+λU_F^2+λV_F^2≤ (1+ϵ)^4



We have with probability 0.99,

    AX̅-B_F^2
    ≤ (1+ϵ)AX-B_F^2
    =0

where the first step follows from Theorem <ref> and the second step follows from the assumption.
Recall that C=UV,
so

    [A, AX̅] - [ A,   B] _F^2+λU_F^2+λV_F^2
    
     =   UV - [ A,   B] _F^2+λU_F^2+λV_F^2≤ (1+ϵ)^4

where the last step uses (<ref>).



However,
if C does not have nice structure,
then we need to apply our procedure Split,
which would introduce the additive error δ.
 Overall, by rescaling ϵ,
our main result is summarized as follows.

Given two matrices A ∈^m × n and B ∈^m × d and λ>0, letting 

    = min_U∈^m × n, V∈^n× (n+d) UV - [ A,   B] _F^2+λU_F^2+λV_F^2,

we have that for any ϵ∈ (0,1), there is an algorithm 
that runs in 

    Õ((A) + (B) + d ·(n/ϵ))
 time and outputs a matrix X ∈^n × d such that there is a matrix A∈^m × n, U∈^m × n and V∈^n× (n+d) satisfying that [A,  AX]-UV_F^2≤δ and

    [A,  AX] - [A,   B] _F+λU_F^2+λV_F^2 ≤ (1+ϵ)  + δ








§ TOY EXAMPLE



We first run our FTLS algorithm on the following toy example, for which we have an analytical solution exactly.
Let A∈ℝ^m× n be A_ii=1 for i=1,⋯,n and 0 everywhere else.
Let B∈ℝ^m× 1 be B_n+1=3 and 0 everywhere else.

The cost of LS is 9,
since AX can only have non-zero entries on the first n coordinates,
so the (n+1)-th coordinate of AX-B must have absolute value 3. Hence the cost is at least 9.
Moreover,
a cost 9 can be achieved by setting X=0 and Δ B=-B.

However,
for the TLS algorithm,
the cost is only 1.
Consider Δ A∈ℝ^m× n where A_11=-1 and 0 everywhere else.
Then C':=[(A+Δ A),   B] does have rank n,
and C'-C_F=1.


For a concrete example,
we set m=10, n=5.
That is, 

    C:=[A,   B]=[ 1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 3; 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0 ]

We first run experiments on this small matrix.
Because we know the solution of LS and TLS exactly in this case,
it is convenient for us to compare their results with that of the FTLS algorithm.
When we run the FTLS algorithm,
we sample 6 rows in all the sketching algorithms.



The experimental solution of LS is C_ which is the same as the theoretical solution. The cost is 9. The experimental solution of TLS is C_ which is also the same as the theoretical result.
The cost is 1.

 

    C_=
    [ 1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0 ]
     C_=
    [ 0 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 3; 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 0 0 0 0 ]


FTLS is a randomized algorithm,
so the output varies.
We post several outputs:

    C_=
    [    0    0    0    0    0    0;    0  0.5  0.5    0    0    0;    0    0    0    0    0    0;    0    0    0    1    0    0;    0    0    0    0  0.1 -0.3;    0    0    0    0 -0.9  2.7;    0    0    0    0    0    0;    0    0    0    0    0    0;    0    0    0    0    0    0;    0    0    0    0    0    0 ]
 
This solution has a cost of 4.3.


    Ĉ_=
    [  0.5 -0.5    0    0    0    0; -0.5  0.5    0    0    0    0;    0    0    0    0    0    0;    0    0 0.09 0.09    0 0.27;    0    0    0    0    0    0;    0    0 0.82 0.82    0 2.45;    0    0    0    0    0    0;    0    0    0    0    0    0;    0    0    0    0    0    0;    0    0    0    0    0    0 ]



This solution has a cost of 5.5455.


    C_=
    [  0.5  0.5    0    0    0    0;    0    0    0    0    0    0;    0    0    1    0    0    0;    0    0    0    0    0    0;    0    0    0    0    1    0;    0    0    0 -0.9    0  2.7;    0    0    0    0    0    0;    0    0    0    0    0    0;    0    0    0    0    0    0;    0    0    0    0    0    0 ]

This solution has a cost of 3.4.
























§ MORE EXPERIMENTS




Figure <ref> shows the experimental result described in Section <ref>.
It collects 1000 runs of our FTLS algorithm on 2 small toy examples.
In both figures, the x-axis is the cost of the FTLS algorithm, measured by C'-C_F^2 where C' is the output of our FTLS algorithm;
the y-axix is the frequency of each cost that is grouped in suitable range.











































