An Analysis of Sketched IRLS for Accelerated
Sparse Residual Regression
Daichi Iwata, Michael Waechter, Wen-Yan Lin, and Yasuyuki Matsushita
Graduate School of Information Science and Technology, Osaka University{iwata.daichi,waechter.michael,lin.daniel,yasumat}@ist.osaka-u.ac.jp

Abstract. This paper studies the problem of sparse residual regression, i.e., learning a linear model using a norm that favors solutions inwhich the residuals are sparsely distributed. This is a common problem in a wide range of computer vision applications where a linear
system has a lot more equations than unknowns and we wish to findthe maximum feasible set of equations by discarding unreliable ones.
We show that one of the most popular solution methods, iteratively
reweighted least squares (IRLS), can be significantly accelerated by theuse of matrix sketching. We analyze the convergence behavior of the
proposed method and show its efficiency on a range of computer vision applications. The source code for this project can be found athttps://github.com/Diwata0909/Sketched_IRLS.
Keywords: sparse residual regression, `1 minimization, randomized algorithm, matrix sketching1

Introduction

We consider the problem of residual minimization, where we wish to learn a linear
model that minimizes the residuals that deviate from the model in some distancemetric. For a linear model we have a matrix A ∈ Rn×d (we consider tall matrices, i.e., the strongly over-determined case with n  d ), a vector b ∈ Rn×1 ,b 6∈ R(A) (the range of A), and we seek to findpx∗ = argmin Ax − b p(1)xfor some p-norm. In this paper, we consider this linear model and call it an `p(residual) minimization problem. For the more general case including non-linear
residual minimization we refer the reader to, e.g., Aftab and Hartley [1] or Kiani
and Drummond [36].
One of the most popular methods for `p residual minimization with 1 ≤ p < 2
is iteratively reweighted least squares (IRLS), in which weighted `2 minimizationis repeatedly computed, eventually converging to the `p solution. It is generally
understood that p = 2 is optimal for Gaussian distributed errors whereas the
1-norm leads to solutions with sparser residuals (or Laplacian distributed errors). For sparse residual regression, the `0 pseudo-norm is appropriate to consider; however, due to its computational complexity, a convex `1 relaxation is2

D. Iwata, M. Waechter, W. Lin, and Y. Matsushita

often employed. Because of its capability of ignoring large outliers and reachingapproximate solutions that yield sparse residuals, `1 residual minimization has
become an important tool for various computer vision tasks. It can be speculated that the core reason of why `2 is still widely preferred over `1 despite `1 ’s
better suitability for many applications, is simply that `2 residual minimization
can be solved efficiently in closed form while `1 residual minimization cannot.
Therefore, acceleration of `1 residual minimization is wanted, especially when
computing budget is limited, e.g., in end-user applications on mobile devices.
In recent years, randomized algorithms are gaining attention beSketched IRLS
cause they offer immense speed up Conventional IRLS
potential with small failure probability in various applications. In
Sketchlinear algebra, matrix operations
can be made efficient with randomized matrix sketching [39,31,58]. It
is based on the idea of approximating an input matrix by multiplying
𝐒(𝐖𝐀)𝐱 ≈ 𝐒(𝐖𝐛)𝐖𝐀𝐱 ≈ 𝐖𝐛a randomly generated sketching
matrix to obtain a much smaller
matrix that still preserves impor- Fig. 1: Illustration of sketched IRLS: Unliketant properties of the input ma- canonical IRLS, we suggest performing antrix. In particular, certain sketch- approximate computation in sketched lower
ing matrices S fulfill the subspace dimensions, yielding significant speed-up
while retaining accuracy.
embedding property [11, Sec. 1.2]
1
2
2
2
Ax 2 ≤ SAx 2 ≤ γ Ax 2γ

(2)

for some γ > 1 with high probability, meaning that the `2 subspace embedding S
preserves the lengths of all vectors Ax within the bounds specified by γ, indicating that we can preserve A’s range via sketching even though the sketchedmatrix SA lives in lower dimensions. This allows us to accelerate `2 regression.
In this paper, we show that `1 residual regression with IRLS can be significantly accelerated by matrix sketching, making it much more useful in practicalapplications. We call the new method sketched IRLS. The key idea is to speed
up the internal computation block by projecting the linear system to lower dimensions using matrix sketching as illustrated in Fig. 1. Through comparisonswith other robust techniques such as RANSAC, we also show that our method
is versatile. This paper’s contributions are summarized as follows:– We propose accelerating IRLS for `1 minimization with matrix sketching,– we analyze the error bound of sketched IRLS compared to canonical IRLS,– and we provide an analysis of our proposed method’s effectiveness on common computer vision problems using both synthetic and real-world data.
The proposed method yields the important benefit that the fundamental computer vision task of regression can be executed in an outlier-robust manner with

An Analysis of Sketched IRLS for Accelerated Sparse Residual Regression

3

the flexible tool of IRLS without excessive computational burden, allowing it to
be used in diverse applications with large datasets.

2

Related works

The problem we study in this paper is related to sparse regression and matrix
sketching. Here we briefly review these two subject areas.
Sparse regression The field of sparse regression has recently developed rapidly,
and its usefulness due to its outlier and noise robustness has been recognized inmany tasks such as face recognition [59,38], image denoising [23], general signal
processing [40], and background estimation [20]. Elad [22] and Zhang [61] givetechnical overviews on sparse regression. The goal of sparse regression is to obtainsolutions where many of the unknowns become zero. A typical setting aims at
minimizing the `0 norm of the solution with an under-determined constraint as
min kxk0
xs.t. Ax = b.

(3)

Similar to this problem is the problem of minimizing the `0 norm of residuals(called robust sensing in [35]) with over-determined Ax ' b, expressed as
min krk0
rs.t.

r = b − Ax.

(4)

Both of these look for sparse unknowns via minimization of the `0 -norm of avector, which is NP-hard. To address this issue, various optimization strategies have been developed. One of them is based on greedy pursuit and tries tofind an approximate solution for `0 minimization directly in an iterative manner [41,47,54,15,43]. Greedy methods are simple and easy to implement, but they
are only guaranteed to work under very restricted conditions [52].
On the other hand, relaxation methods cast the NP-hard `0 minimization
into, for example, `1 minimization, which is a convex surrogate for `0 minimization. Gorodnitsky and Rao proposed the FOCUSS algorithm [28], in which
they approximate `1 minimization with the 2-norm using IRLS. It is also understood that `1 minimization can be solved by linear programming (LP) [26].
Several methods have been proposed to solve the `1 minimization problem, forexample, one based on the interior point method [33], Least Angle Regression(LARS) [21] which is a homotopy algorithm, and the Iterative Shrinkage Thresholding Algorithm (ISTA) [13], to name a few. An excellent summary of recent`1 minimization methods is given by Yang et al . [60].
Among the `1 minimization methods, IRLS is simple and easily implemented
as it only iterates weighted `2 minimization, and it is known for its fast convergence, i.e., only requiring a small number of iterations to achieve an accurate
solution [46,14,9]. Despite the simplicity of the algorithm, `1 minimization with
IRLS shows good performance in recent applications [20,38,42]. Also, IRLS canbe naturally extended to more general `p minimization and minimization withother robust functions such as Huber and the Pseudo-Huber functions [1]. To

4

D. Iwata, M. Waechter, W. Lin, and Y. Matsushita

handle nonlinear problems, generalized IRLS has also been proposed and analyzed theoretically [49,44]. Due to this simplicity and versatility, IRLS has been
one of the most popular methods for `1 minimization and related problems.
Matrix sketching Randomized algorithms are attracting attention recently
as methods to speed up fundamental computation. In linear algebra, matrix
sketching is a randomized algorithm for matrix operations. Woodruff [58], Mahoney [39], Halko et al . [31], and Drineas et al . [16] give technical overviews
on matrix sketching and randomized algorithms. In recent years, randomized
algorithms have been widely applied to some specialized linear and algebraic
problems, e.g., consistent linear systems [29], symmetric diagonally dominantlinear systems [12], and matrix inversion [30]. Matrix sketching is also used foraccelerating several types of matrix decomposition, e.g., singular value decomposition [31], QR decomposition [19], and dynamic mode decomposition [24].
As mentioned, matrix sketching is based on the idea that the range of an
input matrix can be well approximated with high probability by random projection. Such a sketched matrix can then for example be used in regression
problems as they fulfill the subspace embedding property of Eq. (2). One of
the earliest sketching matrices (projectors) were random Gaussian matrices [27]:Let S ∈ Rs×n be a matrix randomly sampled from a Gaussian distribution. Suchan S fulfills the subspace embedding property with γ = 1 + ε and a small ε ≥ 0.
This is a consequence of the Johnson-Lindenstrauss lemma [32]. However, witha dense S sketching takes O(nds) time which is very costly. To overcome this,various methods with other sketching matrices have been proposed. They canroughly be divided into sampling-based and projection-based methods.
Sampling-based methods extract rows of the matrix. The simple intuition isthat, in strongly overdetermined systems, the row vectors are mostly linearly dependent and a subset of them is sufficient to maintain their span. Selecting rowsof the matrix is equivalent to subsampling the data points. Row sampling could
be achieved by premultiplying the input matrix with a binary matrix that picks
the desired rows, but for efficiency this is in practice implemented with simply
selecting those rows in a streaming fashion. The simplest sampling method is
uniform sampling that selects rows with uniform probability. Drineas et al . [17]devised leverage score sampling, which samples based on precomputed leveragescores that unfortunately require computing the singular value decomposition(SVD), making leverage score sampling somewhat impractical.
Projection-based methods premultiply the input matrix with more generalmatrices so that the codomain bases are linear combinations of multiple ofA’s domain bases and not simply selections of A’s domain bases. Ailon andChazelle [3] proposed the Fast Johnson-Lindenstrauss Transform (FJLT), andTropp [53] and Drineas et al . [18] proposed the Subsampled Randomized Hadamard Transform (SRHT) as an extension, which take O(nd log s) time. Clarkson
and Woodruff further proposed the much faster CountSketch method [11], which
was originally used in data streaming [8]. CountSketch takes only O(nnz(A))
time, where nnz(A) is the number of non-zero entries of A.

An Analysis of Sketched IRLS for Accelerated Sparse Residual Regression

35

Proposed method: sketched IRLS

For `p minimization (1 ≤ p < 2) with linear models, IRLS converges to Eq. (1)’ssolution by iteratively minimizing2x(t+1) = argmin W(t) Ay − W(t) b 2(5)ywith a diagonal weight matrix W(t) at the (t)th iteration that is initialized with
W(0) = In . This is called the `p Weiszfeld algorithm [56,57,2]. For `1 minimization, where p = 1, the weights are updated as(t)

Wi,i = Ai,∗ x(t) − bi

− 21

,

where Ai,∗ is matrix A’s ith row, and Wi,i is weight matrix W’s ith diagonal
element. Solving Eq. (5) requires O(nd2 + d3 ) arithmetic operations and is thus
expensive for large matrices. Further, since it is repeatedly solved with updated
W in IRLS, it is the computational bottleneck for IRLS-based `p minimization.
The key idea of the proposed method is to accelerate this computation block,the weighted `2 minimization of Eq. (5), with matrix sketching. Via matrix
sketching, we reduce n-dimensional to much smaller s-dimensional vectors sothat the computational complexity is significantly reduced. Specifically, Eq. (5)’sweighted `2 minimization is modified as2

min WAx − Wb 2xsketch−−−−→g 2,
g − Wbmin WAx2xg ∈ Rs×d , s  n, with retaining the solution’s accuwhere WA ∈ Rn×d and WA
racy. For adopting matrix sketching in IRLS, there are two aspects to consider;(1) when to sketch in the algorithm, and (2) the choice of the sketching method.
3.1

When to sketch?

There are two possible points in time for applying matrix sketching in IRLS; (1)sketch only once before all IRLS iterations (named sketch-once), and (2) sketch
in every iteration of the weighted `2 minimization (named sketch-iteratively).
We analyze the behavior of these two strategies in this paper.
Sketch-once Not only the `2 minimization but also the sketching contributes
to the algorithm’s overall runtime. One way to reduce the sketching costs is to
only perform it once at the beginning of IRLS, outside of IRLS’s iteration loop.
Algorithm 1 shows pseudo code for IRLS with sketch-once.
Sketch-iteratively Unlike sketch-once, sketch-iteratively performs sketching
within the IRLS iteration loop. In each IRLS iteration, we generate an S ∈ Rs×n
with s  n and instead of Eq. (5) we solve(t)(t)
x(t+1) = argmin SW| {z A} x − SW
| {z b}
x
(t)

g
=WA

(t)

g
=Wb

2
.
2

6

D. Iwata, M. Waechter, W. Lin, and Y. Matsushita

Algorithm 1 IRLS with sketch-once
Input: A ∈ Rn×d , b ∈ Rn , sketching size s, #iterations T , threshold δ for termination
Output: Approximate solution x
f (0) ← Identity matrix IsWx(0) ← 0e ← sketch(A, b, s)e bA,for t = 0 : T doef (t) A,e W
f (t) b)x(t+1) ← solve linear least squares(W
(t+1)
(t)
if x
−x
<
δ
then
break
end
if
2
− 1
f (t+1) ← max{, |Ae i,∗ x(t+1) − eWbi |} 2i,i
end for
Algorithm 2 IRLS with sketch-iteratively
Input: A ∈ Rn×d , b ∈ Rn , sketching size s, #iterations T , threshold δ for termination
Output: Approximate solution x
W(0) ← Identity matrix In
x(0) ← 0
for t = 0 : T dog (t) ← sketch(W(t) A, W(t) b, s)
g (t) , Wb
WAg (t) )
g (t) , Wbx(t+1) ← solve linear least squares(WAif x(t+1) − x(t) 2 < δ then break end if
− 1
(t+1)
← max{, |Ai,∗ x(t+1) − bi |} 2
Wi,iend for

While sketch-iteratively requires more computation time than sketch-once because of the sketching operation in the loop, it is expected to be more stablebecause the linear system is sketched differently in each iteration. Pseudo codefor IRLS with sketch-iteratively is shown in Algorithm 2.
3.2

Sketching method choice

In Sec. 2, we discussed various matrix sketching methods. In this study, we
mainly focus on uniform sampling and CountSketch, because of their computational efficiency. Both make the sketching matrix sparse, resulting in sketchingcomputation costs of only O(nnz(A)) time. In practice, they can be implemented
in a streaming fashion without explicitly forming the sketching matrices S.
Uniform sampling samples rows of a tall matrix [A|b] or W[A|b] with uniform probability so that the row-dimension can be reduced. Uniform sampling ofa massively over-determined system has been employed in the past in practicalapplications, in a similar manner to the sketch-once approach described in theprevious section. In an explicit matrix form, the sketching matrix S is a sparse
matrix in which each row has only one element that is 1 while the rest are zeros.
CountSketch [11] is similar to uniform sampling, but instead of subsamplingrows of matrix A, it uses a sketching matrix S in which each column has a single
randomly chosen non-zero entry (typically, 1 or −1).

An Analysis of Sketched IRLS for Accelerated Sparse Residual Regression

4

7

Theoretical analysis

We now analyze the errors in IRLS solutions obtained with and without sketching. To this end, we will show that the proposed sketched IRLS can reliably derivea solution close to the non-sketched case. As mentioned, sketching can approximate a matrix’s range with a certain accuracy. Since “sketch-iteratively” uses
sketching repeatedly, we consider the errors introduced in each iteration. The
goal is to reveal the relationship between Aex − b 1 and Ax∗ − b 1 for IRLS,∗
e is the approximate solution obtained
where x is the optimal solution, and x
with sketched IRLS. We derive it by combining the error bounds of sketchinge(t) be the solutions of canonical andand IRLS’s convergence rate. Let x∗(t) and xe∗(t+1) be the solution obtained bysketched IRLS, resp., after t iterations, and xe(t) .
solving Eq. (5) without sketching and using W based on x
4.1

Condition for the residual to decrease

Before considering error bounds, we first show what condition must be fulfilledx(t+1) −so that the residual in sketched IRLS decreases monotonically, i.e., Ae(t)x − b 1 . For matrix sketching with `2 regression problems, severalb 1 ≤ Aeerror bounds are known [39,48,17]. In a general form, we can write them asAex − b 2 ≤ (1 + ε) Ax∗ − b 2 ,
e−x
x

∗
2

= ∆x 2 ≤ εx ,

(6)
(7)

e − x∗ := ∆x and small errors ε and εx , which depend on the sketching
using xmethod and sketching size. The error decay rate is known for canonical IRLS,
e.g., linear and super-linear convergence rates [49,14]. According to Daubechies [14,Sec. 6.1], for sparse variable x problems (Eq. (3)) we havee(t) − x∗ 1e∗(t+1) − x∗ 1 ≤ µ xx(8)with a constant µ ≤ 1. It is known that sparse variable problems (Eq. (3)) andsparse residual problems (Eq. (4)) can be treated as the same problem underx = r [7], and for sparse residual problems we have(Aex∗(t+1) − b) − (Ax∗ − b) 1 ≤ µ (Aex(t) − b) − (Ax∗ − b) 1 .

(9)

From these, the residuals Aex(t+1) − b 1 and Aex(t) − b 1 satisfy(Aex(t+1) − b) 1 − (Ax∗ − b) 1 ≤ (Aex(t+1) − b) − (Ax∗ − b) 1
= (Aex∗(t+1) − b) + A∆x − (Ax∗ − b) 1
≤ (Aex∗(t+1) − b) − (Ax∗ − b) 1 + A∆x 1
Eq. (9)≤ µ (Aex(t) − b) − (Ax∗ − b) 1 + A∆x 1≤ µ (Aex(t) − b) 1 + µ (Ax∗ − b) 1 + A∆x 1 .

(10)

8

D. Iwata, M. Waechter, W. Lin, and Y. Matsushita

From Eq. (10), we finally have(Aex(t+1) − b) 1 ≤ µ (Aex(t) − b) 1 + (1 + µ) (Ax∗ − b) 1 + A∆x 1 ,(11)and know that when A∆x 1 ≤ (1 − µ) (Aex(t) − b) 1 − (1 + µ) (Ax∗ − b) 1
holds, then Ae
x(t+1) −b 1 ≤ Ae
x(t) −b 1 . We remark that ∆x satisfies Eq. (7),and when the error of sketching ∆x is sufficiently close to 0, then sketched IRLS
lets the objective decrease monotonically.
4.2

Worst-case residual boundIn this section, we show the relationship between the residuals of canonical IRLS
x(t) − b 1 after t IRLS iterations using
Ax∗(t) − b 1 and sketched IRLS Aenorm relations. For canonical IRLS, after t iterations we have
η Ax∗(0) − b 1 = Ax∗(t) − b 1 ,(12)where, assuming that t is big enough for x∗(t) to converge to x∗ , η ∈ [0, 1] is theratio between the `1 residuals of the solutions for Eq. (1) with p = 1 and p = 2.
For sketched IRLS, under Section 4.1’s condition, after t iterations we havex(0) − b 1 .
Aex(t) − b 1 ≤ Ae

(13)

Further, we have the norm relations
kxk2 ≤ kxk1 ≤√n kxk2 ,(14)

where n is the number of dimensions of x. We can now derive a bound on thesolution error due to the sketching after t IRLS iterations:Eq. (14)√
Eq. (6)√
Eq. (13)Ae
x(t) − b 1 ≤ Ae
x(0) − b 1 ≤ n Aex(0) − b 2 ≤ n (1+ε) Ax∗(0) − b 2√

Eq. (14)≤

Eq. (12)n (1+ε) Ax∗(0) −b 1 =√n (1+ε)η −1 Ax∗(t) −b 1 . (15)This bound may not be very tight if n is large; however, we next show that atighter bound can be expected in practice.
4.3

Expected residual boundThe reason why the bound shown in Eq. (15) looks somewhat loose is mainly
due to the right side expression of Eq. (14). The bound expresses the worst caseand if that case rarely occurs, the bound does not have strong implications.
Therefore, it is reasonable to consider the expected value for the `1 -norm. Here,let x ∈ Rn be an arbitrary vector with kxk2 = r and the expected value forthe `1 -norm be E[kxk1 ]. The expectation of the `1 -norm of a vector x becomesE[kxk1 ] = E[|x1 |] + . . . + E[|xn |]. In polar coordinates, x1 , . . . , xn are


x1




 x2

= r cos θ1 ,= r sin θ1 cos θ2 ,···xn−1 = r sin θ1 sin θ2 . . . sin θn−2 cos θn−1 ,



 xn
= r sin θ1 sin θ2 . . . sin θn−2 sin θn−1 .

An Analysis of Sketched IRLS for Accelerated Sparse Residual Regression

9Here, we assume that arbitrary random vectors are generated by uniformly
distributed θ, which means that a probability density function p(θ) is a constant. In this case, the probability density function for x1 depends on the arcsine distribution. Namely, for −r ≤ x1 ≤ r, the probability density function isp(x1 ) = √ 12 2 . Therefore, the expectation E[|x1 |] becomes
π

r −x1

Z r
E[|x1 |] =Z r|x1 | p(x1 ) dx1 =−r1
2|x1 | pdx1 =2
ππ r 2 − x1
−r

Z r
0

x1
p

r2 − x21

 q
r
2
2
2
2
− r − x1 = r.
=
π
π
0

dx1
(16)

We can write x2 as x2 = r1 cos θ2 with r1 = r sin θ1 , and we can obtain E[|x2 |] by
R E[r ]
using Eq. (16) recursively, i.e., E[|x2 |] = −E[r1 1 ] |x2 |p(x2 ) dx2 = π2 E[r1 ] = ( π2 )2 r.
Finally, the expected value E[kxk1 ] becomesE[kxk1 ] = E[|x1 |] + E[|x2 |] + . . . + E[|xn−1 |] + E[|xn |]
2
n−1
n−1
r + π2r= π2 r + π2 r + . . . + π2
n−1
2
n−1
P 2 i+( 2 )n−1 −2( 2 )nr + π2r = π π 1− 2 π r.
=
π
π

i=1

We thus have an expected bound of
−1

∗(t)

Ae
x(t) − b 1 ≤

2 n−1
2 n
2
−2( π
)
π +( π )
2
1− π

(1 +

2
'
. With n → ∞, the fraction expression converges to π−2
1
√

ε) η
Ax
−b
1.75. The expected value is considerably smaller than the worst case value n,
indicating that when n is large, Sec. 4.2’s worst case bound is hardly relevant.
Although it is still far from a strict bound due to the dependency on η, as we
will see in the following sections, sketched IRLS yields highly accurate approximations in practice in various settings.

5

Performance evaluation

We first evaluate the proposed method’s performance in comparison to canonical
IRLS using synthetic data in common problems, namely, in residual minimizationand low-rank approximation.
5.1

Residual minimization

In `1 residual minimization, given A ∈ Rn×d , n > d and b ∈ Rn , we wish to solve
min kAx − bk1x

for x ∈ Rd . To form a highly over-determined linear system, we set the size ofmatrix A to (n, d) = (106 , 40). For assessing the performance variations with respect to the outlier distribution, we formed matrix A in three distinct ways: (a) A10

D. Iwata, M. Waechter, W. Lin, and Y. Matsushita
uniform sampling, Sketch-Once
uniform sampling, Sketch-Iteratively

leverage score sampling, Sketch-Onceleverage score sampling, Sketch-Iteratively

CountSketch, Sketch-Once
CountSketch, Sketch-Iteratively

SRHT, Sketch-Once
SRHT, Sketch-Iteratively

no sketching
RANSAC



HUURU

HUURU

HUURU















WLPHVHF





(a) Uniform data 20 %




















WLPHVHF





(b) Uniform data 60 %













WLPHVHF









(c) Biased data

Fig. 2: Averages of the error d1 kx∗ − x(t) k2 and the standard deviations over timein residual minimization with various matrix sketching methods and RANSAC
on uniform synthetic data with a lower outlier rate (left), with a higher outlier
rate (center ), and with biased synthetic data (right).

drawn from a uniform distribution in [0, 10] and flipped signs of 20 % elements tocreate outliers (called uniform data 20 %, hereafter), (b) A created like (a) but
with 60 % outliers (called uniform data 60 %, hereafter), and (c) A created like(a) but further corrupted by adding a large value (= 103 ) to randomly selected
0.1 % of rows (called biased data, hereafter). The ground truth x∗ is created,and based on x∗ , b is pre-computed before adding outliers. We evaluate theerror defined as d1 kx∗ − x(t) k2 , where x(t) is the solution after the tth iteration.
We also compare the accuracy with RANSAC [25] to assess the robustness of `1
minimization against outliers. Unlike `1 minimization, RANSAC requires adjusting a few parameters. We changed the number of samplings, RANSAC’s mostimportant parameter, and chose the best result from {d + 1, d × 10, d × 102 }.
Figures 2a, 2b and 2c show the results of uniform data 20 %, uniform data
60 %, and biased data, resp. All results are averages of 10 trials with different
random seeds. The plots show averages of the error, error bars indicate standard
deviations. From Figs. 2a and 2b we can observe that `1 minimization works wellfor problems with relatively few outliers, and RANSAC shows slow convergences.
Both methods need to solve least squares minimization many times, IRLS canbe expected to improve the solution for each loop, but RANSAC does not necessarily do so. RANSAC is further known to fail at high-dimensional problems.
Regarding sketching methods, while sampling-based methods with sketch-oncework well for uniform data as shown in Fig. 2a, their accuracies become lower andvariances become larger on the biased data as shown in Fig. 2c. Projection-based
methods work well for either data with low variances, especially CountSketchwith sketch-iteratively. For leverage score sampling, the comparison was doneusing the known leverage scores while it is usually unknown a priori. From hereon out, we will focus only on uniform sampling and CountSketch as the chosen
sketching methods because of their efficiency as discussed in Sec. 3.2.

An Analysis of Sketched IRLS for Accelerated Sparse Residual Regression
CountSketch, Sketch-Once
CountSketch, Sketch-Iteratively
















WLPHVHF

H

UDQN





ALM
ALM with FRSVTH 



HUURU

HUURUH 

no sketching   uniform sampling, Sketch-Once
uniform sampling, Sketch-Iteratively   11





VNHWFKLQJVL]HUDQN





WLPHVHF


H

Fig. 4: Comparison between
sketch-iteratively and sketchFig. 3: Error over time in `1 singular value de- once. Colors indicate fasterfor
sketchcomposition on synthetic data with lower outlier convergence
iterativelyorsketch-once.
rate, (left), with higher outlier rate (right)(a) 10 % outliers

5.2(b) 40 % outliers

Low-rank approximation

Next, we evaluate the proposed method on low-rank approximation with `1 singular value decomposition. Given a matrix M ∈ Rm×n , its rank-r approximation(r < min(m, n)) with `1 -norm can be written as
min kM − UV> k1 ,U,Vwhere U and V> are m×r and r ×n matrices, respectively. The solution method
involves `1 minimization subproblems that are computed iteratively [34] asU ← argmin kM − UV> k1 ,U

V ← argmin kM − UV> k1 ,Vstarting from a random initialization of V.
44
We generated an M ∈ R10 ×10 with rank(M) = 40, flipped signs of 10 % and
40 % of the elements. We set the sketching size s as 1, 600 and 3, 200. We also
compare the accuracy with robust principal component analysis (R-PCA) [6].
R-PCA decomposes a corrupted matrix M into a low-rank matrix A and asparse corruption matrix E, i.e., M = A + E. To solve R-PCA, the Augmented
Lagrange Multiplier method (ALM) [37] is a known effective approach. Also aprevious study accelerated ALM by fast randomized singular value thresholding(FRSVT) [45]. We used native ALM and ALM with FRSVT as comparison.
These methods also require tuning hyper-parameter λ, and we chose the best
result from λ ∈ {10−1 , 10−2 , 10−3 }.
1kM∗ − M(t) kF , where M∗
To assess the accuracy, we define the error as mn(t)
is the uncorrupted original matrix and M is the estimated low-rank matrix
after the tth iteration. Figures 3a and 3b show the results. In this setting, `1
minimization achieves high accuracy in both datasets. R-PCA converges quickly
but does not work well for the dataset with many outliers. The sketch-once
strategy shows about 3 times faster convergence compared to canonical IRLS in
Fig. 3a, and also uniform sampling strategy shows fast convergence in Fig. 3b.

12

D. Iwata, M. Waechter, W. Lin, and Y. Matsushita



left imageℓ1 min. with Sketched IRLSUHVLGXDOXQLIRUPVDPSOLQJ6NHWFK2QFH
XQLIRUPVDPSOLQJ6NHWFK,WHUDWLYHO\
&RXQW6NHWFK6NHWFK2QFH
&RXQW6NHWFK6NHWFK,WHUDWLYHO\QRVNHWFKLQJ







right imageℓ2 minimization(a) Left: Input images. Right: `1 and `2
stitching results.





WLPHVHF

(b) `1 residual variation over time

Fig. 5: Performance of image stitching with `1 and `2 homography estimation

5.3

When to sketch?

In the following, we show an experiment that serves to give a first intuition
for when a user should pick the sketch-once or the sketch-iteratively sketching
44
regime. We generated R10 ×10 matrices with 10 % noise, we picked the ranks
from {10, 20, 40, 80} and sketching sizes from {rank×5, ×10, ×20, ×40, ×80}, andconducted a low-rank approximation experiment. We checked the times tSO andtSI until sketch-once and sketch-iteratively achieved an accuracy of 10−5 . Figure 4 shows the values of ln(tSO /tSI ). For each rank, the larger the sketching size,
the faster sketch-once converges. As the matrix rank increases, sketch-iteratively
shows faster convergence at larger sketching sizes. The sketching size for ideal
approximations increases faster than O(r). When the sketching size is not big
enough, the risk of not making ideal approximations becomes higher. Especially,at small sketching sizes, sketch-once is susceptible to making bad approximations, whereas sketch-iteratively alleviates this effect by repeated sketching andshows faster convergence for smaller sketching sizes.

6

Applications

Countless computer vision tasks are built upon robust regression. In practicalscenarios, the input is typically quite large and noisy, either due to a noisy
capturing process or a noisy precomputation step, e.g., incorrect feature matches.
In this section, we demonstrate the proposed sketched IRLS’s effectiveness onsuch real-world applications. We adapted our method to problems with overdetermined system: homography estimation and point cloud alignment.
6.1

Homography estimation

Homography estimation is an important building block for image stitching andplane-to-plane alignment. Given correspondences in two images, the image trans-

An Analysis of Sketched IRLS for Accelerated Sparse Residual Regression

13

formation can be written by a 3×3 homography H with 8 DoF as
λ[x0 , y 0 , 1]> = H[x, y, 1]> ,where λ is a scaling factor, and (x, y) and (x0 , y 0 ) are corresponding points in theleft and right images, respectively [50]. Given a set of correspondences, one canestimate H by solving a homogeneous system [51, Chap. 6.1]. The set of correspondences may contain many wrong correspondences from erroneous feature
matching. Therefore, conventional methods use robust estimation techniques,
such as RANSAC or `1 regression.
Result We obtained point correspondences by matching AKAZE [4] featuresin the two images shown in the left column of Fig. 5a. From 17, 076 obtained
point correspondences we estimated the homography. The second column ofFig. 5a shows the stitching results of `1 and `2 minimization, respectively. `1minimization successfully joined the two images whereas `2 minimization produced a strongly erroneous result due to feature mismatches. In this experiment,
we sketched the matrix to 17, 076×2×1 % ' 341 equations (note that each point
pair gives two equations), and confirmed that the solution did not differ fromthe direct `1 minimization without sketching. In Fig. 5b, we can see that the
proposed method (uniform sampling + sketch-once) converges about 5× fasterthan canonical IRLS while maintaining the same accuracy.
6.2

Point cloud alignment

Consider the two point sets P = [p1 , . . . , pl ] and Q = [q1 , . . . , ql ] captured fromdifferent viewpoints around an object. For each i, the points pi ∈ R3 and qi ∈ R3
approximately correspond to the same surface point on the 3D object. Since
P may be positioned and scaled differently in space than Q, we search for asimilarity transform T that, applied to all points in P, makes P and Q roughly
coincide. We therefore optimize
X
min
Tp̃i − q̃i ,
(17)
T

i

with the tilde denoting homogeneous representations. This problem is commonly
solved with `2 minimization but with large outliers, e.g., due to wrong point
correspondences, `1 minimization may be superior. In the veteran but still frequently used Iterative Closest Point (ICP) algorithm [10,5], the point sets are
first manually pre-aligned, then it searches for the nearest neighbor in Q for eachpoint in P, optimize Eq. (17), and iterate.
Result In this task, we used the Stanford Bunny [55]. The two input point
clouds are shown in Fig. 6a (left). Each set consists of 105 points, viewed fromdifferent viewpoints in different scalings. We applied the proposed sketched IRLS
to perform ICP with `1 residual minimization. The second and third figures ofFig. 6a show the results of `1 minimization with sketched IRLS and conventional`2 minimization. It is observed that `1 minimization results in an accurate alignment of the two point clouds, unaffected by inaccurate correspondences.

14

D. Iwata, M. Waechter, W. Lin, and Y. Matsushita

HUURU

XQLIRUPVDPSOLQJ6NHWFK2QFH
XQLIRUPVDPSOLQJ6NHWFK,WHUDWLYHO\
&RXQW6NHWFK6NHWFK2QFH
&RXQW6NHWFK6NHWFK,WHUDWLYHO\QRVNHWFKLQJ

initial poseℓ1 min. with
Sketched IRLS

ℓ2 minimization(a) Input point clouds and alignment results

10 3









WLPHVHF



(b) Error variation over time

Fig. 6: Performance of point cloud alignment via ICP with `1 and `2 similarity
transform estimation

For this experiment we set the sketching size s as 10 % of the original problem,and gradually transformed the point set Q with fixed P. Since we have theground truth here, we evaluate the error by 1l kQ∗ − Q(t) kF , where Q(t) is theresult after the tth iteration. The evaluation of the error variation w.r.t. time
only for computing the transformation (excluding matchings) is summarized
in Fig. 6b. While CountSketch + sketch-once did not show good convergence,
the other sketching methods find a good alignment with significant speed upcompared to the conventional method.

7

Discussion

Our experiments showed that sketching is effective in reducing IRLS’s runtime
and `1 minimization with IRLS works well on a wide range of computer vision
problems. Other robust methods such as RANSAC and RPCA are certainly goodat specific problems, but IRLS `1 minimization is versatile without requiring
parameter tuning and its convergence is demonstrably superior in some tasks.
Regarding when to sketch, sketch-once is superior if the rank r of the design
matrix A is very small and the sketching is not aggressive, i.e., s  r. However,if the rank is high or we sketch aggressively, e.g., s < 5r, then it is likely thatA’s range will not be preserved, and we need to perform sketch-iteratively to beable to recover from badly chosen samples.
If one naı̈vely performs subsampling on the input data, this would be equal
to sketch-once sketching, basically treating IRLS as a black box that can never
be touched again after the data has initially been subsampled. The experiments
showed that in applications where sketch-iteratively performs better, we want
to open that black box and perform subsampling in every iteration.
Acknowledgments: This work is supported by JSPS CREST Grant Number
JPMJCR1764, Japan. Michael Waechter was supported through a postdoctoral
fellowship by the Japan Society for the Promotion of Science (JP17F17350).

An Analysis of Sketched IRLS for Accelerated Sparse Residual Regression

15

References
1. Aftab, K., Hartley, R.: Convergence of iteratively re-weighted least squares to robust M-estimators. In: Winter Conference on Applications of Computer Vision(WACV) (2015) 1, 3
2. Aftab, K., Hartley, R., Trumpf, J.: Generalized Weiszfeld algorithms for Lq optimization. Transactions on Pattern Analysis and Machine Intelligence (PAMI)
37(4), 728–745 (2015) 5
3. Ailon, N., Chazelle, B.: Approximate nearest neighbors and the fast JohnsonLindenstrauss transform. In: Symposium on Theory of Computing (2006) 4
4. Alcantarilla, P.F., Nuevo, J., Bartoli, A.: Fast explicit diffusion for accelerated
features in nonlinear scale spaces. In: British Machine Vision Conference (BMVC)(2013) 13
5. Besl, P.J., McKay, N.D.: A method for registration of 3-D shapes. Transactions on
Pattern Analysis and Machine Intelligence (PAMI) 14(2), 239–256 (1992) 13
6. Candès, E.J., Li, X., Ma, Y., Wright, J.: Robust principal component analysis?
Journal of ACM 58(3), 11:1–11:37 (2011) 11
7. Candes, E.J., Tao, T.: Decoding by linear programming. IEEE Transactions on
Information Theory 51(12), 4203–4215 (2005) 7
8. Charikar, M., Chen, K., Farach-Colton, M.: Finding frequent items in data streams.
In: International Colloquium on Automata, Languages and Programming (2002)
4
9. Chen, C., He, L., Li, H., Huang, J.: Fast iteratively reweighted least squares algorithms for analysis-based sparse reconstruction. Medical Image Analysis 49, 141 –
152 (2018) 3
10. Chen, Y., Medioni, G.: Object modeling by registration of multiple range images.
In: International Conference on Robotics and Automation (ICRA) (1991) 13
11. Clarkson, K.L., Woodruff, D.P.: Low rank approximation and regression in input
sparsity time. In: Symposium on Theory of Computing (2013) 2, 4, 6
12. Cohen, M.B., Kyng, R., Miller, G.L., Pachocki, J.W., Peng, R., Rao, A.B., Xu,
S.C.: Solving SDD linear systems in nearly m log1/2 n time. In: Symposium on
Theory of Computing (2014) 4
13. Daubechies, I., Defrise, M., De Mol, C.: An iterative thresholding algorithm forlinear inverse problems with a sparsity constraint. Communications on Pure andApplied Mathematics 57(11), 1413–1457 (2004) 3
14. Daubechies, I., DeVore, R., Fornasier, M., Güntürk, C.S.: Iteratively reweighted
least squares minimization for sparse recovery. Communications on Pure and Applied Mathematics 63(1), 1–38 (2010) 3, 7
15. Donoho, D.L., Tsaig, Y., Drori, I., Starck, J.L.: Sparse solution of underdetermined
systems of linear equations by stagewise orthogonal matching pursuit. Transactions
on Information Theory 58(2), 1094–1121 (2012) 3
16. Drineas, P., Mahoney, M.W.: Randnla: Randomized numerical linear algebra. Communications of the ACM 59(6), 8090 (2016) 4
17. Drineas, P., Mahoney, M.W., Muthukrishnan, S.: Sampling algorithms for `2 regression and applications. In: ACM-SIAM Symposium on Discrete Algorithm (2006)
4, 7
18. Drineas, P., Mahoney, M.W., Muthukrishnan, S., Sarlós, T.: Faster least squares
approximation. Numerische Mathematik 117(2), 219–249 (2011) 4
19. Duersch, J., Gu, M.: Randomized QR with column pivoting. SIAM Journal on
Scientific Computing 39(4), 263–291 (2017) 4

16

D. Iwata, M. Waechter, W. Lin, and Y. Matsushita

20. Dutta, A., Richtárik, P.: Online and batch supervised background estimation via
l1 regression. In: Winter Conference on Applications of Computer Vision (WACV)(2019) 3
21. Efron, B., Hastie, T., Johnstone, I., Tibshirani, R.: Least angle regression. Annalsof Statistics 32(2), 407–499 (2004) 3
22. Elad, M.: Sparse and Redundant Representations: From Theory to Applicationsin Signal and Image Processing. Springer, 1st edn. (2010) 3
23. Elad, M., Aharon, M.: Image denoising via sparse and redundant representations
over learned dictionaries. Transactions on Image Processing 15(12), 3736–3745(2006) 3
24. Erichson, N., Mathelin, L., Brunton, S., Kutz, J.: Randomized dynamic mode
decomposition. SIAM Journal on Applied Dynamical Systems 18(4), 1867–1891(2017) 4
25. Fischler, M.A., Bolles, R.C.: Random sample consensus: A paradigm for model
fitting with applications to image analysis and automated cartography. Communications of the ACM 24(6), 381–395 (1981) 10
26. Gentle, J.E.: Matrix Algebra: Theory, Computations, and Applications in Statistics. Springer (2007) 3
27. Goldstine, H.H., von Neumann, J.: Numerical inverting of matrices of high order.
ii. In: American Mathematical Society (1951) 4
28. Gorodnitsky, I.F., Rao, B.D.: Sparse signal reconstruction from limited data using FOCUSS: A re-weighted minimum norm algorithm. Transactions on Signal
Processing 45(3), 600–616 (1997) 3
29. Gower, R., Richtárik, P.: Randomized iterative methods for linear systems. SIAM
Journal on Matrix Analysis and Applications 36(4), 1660–1690 (2015) 4
30. Gower, R., Richtrik, P.: Randomized quasi-Newton updates are linearly convergent
matrix inversion algorithms. SIAM Journal on Matrix Analysis and Applications
38(4), 1380–1409 (2016) 4
31. Halko, N., Martinsson, P.G., Tropp, J.A.: Finding structure with randomness:
Probabilistic algorithms for constructing approximate matrix decompositions.
SIAM Rev. 53(2), 217–288 (2011) 2, 4
32. Johnson, W.B., Lindenstrauss, J.: Extensions of Lipschitz mappings into a Hilbert
space. Contemporary Mathematics 26, 189–206 (1984) 4
33. Karmarkar, N.: A new polynomial-time algorithm for linear programming. Combinatorica 4(4), 373–395 (1984) 3
34. Ke, Q., Kanade, T.: Robust L1 norm factorization in the presence of outliers andmissing data by alternative convex programming. In: Conference on Computer
Vision and Pattern Recognition (CVPR) (2005) 11
35. Kekatos, V., Giannakis, G.B.: From sparse signals to sparse residuals for robust
sensing. IEEE Transactions on Signal Processing 59(7), 3355–3368 (2011) 3
36. Kiani, K.A., Drummond, T.: Solving robust regularization problems using iteratively re-weighted least squares. In: Winter Conference on Applications of Computer Vision (WACV) (2017) 1
37. Lin, Z., Chen, M., Wu, L., Ma, Y.: The augmented Lagrange multiplier method forexact recovery of corrupted low-rank matrices. Tech. Rep. UILU-ENG-09-2215, Coordinated Science Laboratory, University of Illinois at Urbana-Champaign (2010)11
38. Lu, C., Lin, Z., Yan, S.: Smoothed low rank and sparse matrix recovery by iteratively reweighted least squares minimization. Transactions on Image Processing
24(2), 646–654 (2015) 3

An Analysis of Sketched IRLS for Accelerated Sparse Residual Regression

17

39. Mahoney, M.W.: Randomized algorithms for matrices and data. Foundations andTrends R in Machine Learning 3(2), 123–224 (2011) 2, 4, 7
40. Mallat, S.G.: A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way.
Academic Press, Inc., 3rd edn. (2008) 3
41. Mallat, S.G., Zhang, Z.: Matching pursuits with time-frequency dictionaries. Transactions on Signal Processing 41(12), 3397–3415 (1993) 3
42. Millikan, B., Dutta, A., Rahnavard, N., Sun, Q., Foroosh, H.: Initialized iterative
reweighted least squares for automatic target recognition. In: Military Communications Conference (2015) 3
43. Needell, D., Tropp, J.A.: CoSaMP: Iterative signal recovery from incomplete andinaccurate samples. Applied and Computational Harmonic Analysis 26(3), 301–
321 (2009) 3
44. Ochs, P., Dosovitskiy, A., Brox, T., Pock, T.: On iteratively reweighted algorithms
for non-smooth non-convex optimization in computer vision. SIAM Journal on
Imaging Sciences 8(1), 331–372 (2015) 4
45. Oh, T.H., Matsushita, Y., Tai, Y.W., Kweon, I.S.: Fast randomized singular value
thresholding for low-rank optimization. Transactions on Pattern Analysis and Machine Intelligence (PAMI) 40(2), 376–391 (2018) 11
46. Osborne, M.R.: Finite Algorithms in Optimization and Data Analysis. John Wiley& Sons, Inc. (1985) 3
47. Pati, Y.C., Rezaiifar, R., Krishnaprasad, P.S.: Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition. In:
Asilomar Conference on Signals, Systems, and Computers (1993) 3
48. Sarlós, T.: Improved approximation algorithms for large matrices via random projections. In: Symposium on Foundations of Computer Science. pp. 143–152 (2006)
7
49. Sigl, J.: Nonlinear residual minimization by iteratively reweighted least squares.
Computational Optimization and Applications 64(3), 755–792 (2016) 4, 7
50. Szeliski, R.: Video mosaics for virtual environments. Computer Graphics and Applications 16(2), 22–30 (1996) 13
51. Szeliski, R.: Computer Vision - Algorithms and Applications. Texts in Computer
Science, Springer (2011) 13
52. Tropp, J.A.: Greed is good: Algorithmic results for sparse approximation. Transactions on Information Theory 50(10), 2231–2242 (2004) 3
53. Tropp, J.A.: Improved analysis of the subsampled randomized Hadamard transform. Advances in Adaptive Data Analysis 3, 115–126 (2011) 4
54. Tropp, J.A., Gilbert, A.C., Strauss, M.J.: Algorithms for simultaneous sparse approximation. Part I. Signal Process. 86(3), 572–588 (2006) 3
55. Turk, G., Levoy, M.: Zippered polygon meshes from range images. In: SIGGRAPH(1994) 13
56. Weiszfeld, E.: Sur le point pour lequel la somme des distances de n points donnésest minimum. Tohoku Mathematical Journal (1937) 5
57. Weiszfeld, E., Plastria, F.: On the point for which the sum of the distances to n
given points is minimum. Annals of Operations Research (2009) 5
58. Woodruff, D.P.: Sketching as a tool for numerical linear algebra. Foundations andTrends R in Theoretical Computer Science 10(1-2), 1–157 (2014) 2, 4
59. Wright, J., Yang, A.Y., Ganesh, A., Sastry, S.S., Ma, Y.: Robust face recognition
via sparse representation. Transactions on Pattern Analysis and Machine Intelligence (PAMI) 31(2), 210–227 (2009) 3

18

D. Iwata, M. Waechter, W. Lin, and Y. Matsushita

60. Yang, A., Zhou, Z., Ganesh Balasubramanian, A., Sastry, S., Ma, Y.: Fast L1minimization algorithms for robust face recognition. Transactions on Image Processing 22(8), 3234–3246 (2013) 3
61. Zhang, Z., Xu, Y., Yang, J., Li, X., Zhang, D.: A survey of sparse representation:
Algorithms and applications. IEEE Access 3, 490–530 (2015) 3