Robust Least Squares for Quantized Data Matrices
Richard Clancy Stephen Becker
math.OC cs.NA math.NA 65K05, 90C99



 


stephen.becker@colorado.edu

[cor1]Corresponding author

richard.clancy@colorado.edu



Department of Applied Mathematics at the University of Colorado, Boulder, CO 80309




	In this paper we formulate and solve a robust least squares problem for a system of linear equations subject to quantization error in the data matrix. Ordinary least squares fails to consider uncertainty in the operator, modeling all noise in the observed signal. Total least squares accounts for uncertainty in the data matrix, but necessarily increases the condition number of the operator compared to ordinary least squares. Tikhonov regularization or ridge regression is frequently employed to combat ill-conditioning, but requires parameter tuning which presents a host of challenges and places strong assumptions on parameter prior distributions. The proposed method also requires selection of a parameter, but it can be chosen in a natural way, e.g., a matrix rounded to the 4th digit uses an uncertainty bounding parameter of 0.5e-4. We show here that our robust method is theoretically appropriate, tractable, and performs favorably against ordinary and total least squares. 



	Least squares approximation Quantization error Ridge regression Robust optimization Subgradient methods Tikhonov regularization Total least squares







§ INTRODUCTION

The primary goal of this paper is to recover unknown parameters from a noisy observation and an uncertain linear operator. In particular, our mathematical model for robust least squares is

    min_x  {max_∈   (+)x -b ^2 }

which we refer to as our robust optimization (RO) problem. The Euclidean norm is denoted by · and 𝒰 is the uncertainty set from which perturbations in A are drawn.  

The above RO formulation is motivated by two situations. In both cases, we let A̅ and x̅ represent the true and unknown data matrix and parameter vector, respectively. We model b = A̅x̅ + η, with η i.i.d. Gaussian, and we only have knowledge of A = A̅ -. We don't know  explicitly but can make inferences based on the problem. In the first situation, we consider a data matrix subject to quantization or round-off error. Suppose the observed matrix A has elements rounded to the hundredth place. Our uncertainty set can be written as = {∈^m × n  : _∞≤δ} with δ = 0.005. If A_i,j = 0.540, then we know the true A̅_i,j∈ (0.535,  0.545], hence _i,j∈ (-0.005, 0.005]. The norm ·_∞ takes the maximum absolute value of any element in the matrix.

The second problem considers a data matrix with uncertainty proportional to the magnitude of each entry, i.e. = {∈^m × n : _i,j∈ (-p |_i,j|, p |_i,j|  ]   }. Here, p denotes a proportionality constant. Data subject to ± 1 % uncertainty would have p = 0.01. 
The two cases cover the effects of finite-precision in fixed and floating point representations, respectively. In both problems, the uncertainty sets are specified element-wise allowing us to decouple along rows. 

Due to limitations in both ordinary (OLS) and total least squares (TLS), the signal processing community has sought alternatives when operator uncertainty is present. In <cit.> and <cit.>, the authors derive maximum likelihood estimators (MLE) and Cramér Rao bounds for problems with data matrices subject to Gaussian noise under two different models; errors-in-variables and a random variable model. The setup for <cit.> only has access to sign measurements which is more restrictive than we consider here. Both papers focus primarily on the case where variance in A and b are known; 
in contrast, we make no distributional assumption on A (it could be stochastic, deterministic, or adversarial). In the absence of knowing the variance, <cit.> shows their estimator reduces to the OLS solution. The problem is treated through an approximate message passing framework in <cit.> for more general, structured perturbations. Our use case fits under this umbrella but they impose a sparsity inducing prior. Although our method performs well with sparse solutions, we do not assume it.

Note that the model in Eq. <ref> differs from other robust least squares problems considered in <cit.>, <cit.> and classic papers from the late 1990's <cit.> and early 2000's <cit.>. Those works usually made special assumptions like the constraint and objective norms match (e.g. minimize x_p subject to Ax- b_p ≤ v), column-wise separability, or ellipsoidal uncertainty sets. The work in <cit.> is similar regarding row-wise separability, but they impose a hard constraint on the 2-norm of the solution, i.e., x≤ k for k ∈. For an overview of robust optimization, see <cit.> and <cit.> with the latter focusing on implementation.

In this paper, we consider box constraints over entries of the data matrix. Our main contribution is the formulation of a robust objective to handle quantization error in least squares problems and the presentation of methods to solve it. Although the proposed method requires selection of an uncertainty parameter, it is chosen in a natural and theoretically appropriate way based on the observed extent of quantization. This is in contrast to ridge regression and MLE based methods that require involved parameter tuning or a priori knowledge of the probability distributions from which model uncertainty is drawn. We anticipate our method to be most effective under moderate to heavy quantization where fidelity loss is greater than 0.1%.



 §.§ Limitations of Ordinary and Total Least Squares
 
The ordinary least squares (OLS) problem seeks a vector x to minimize the residual given by

    min_x  x-b^2.

We focus our attention on the over-determined case where m>n and further assume that A is full rank. It is well known that the closed form solution to (<ref>) is given by 

    x̂_OLS = (^T)^-1A^T b.

A key assumption for OLS is that A is known with certainty and b is subject to additive noise. 
The OLS solution (<ref>) is the MLE for the model Ax-b = η∼𝒩(0, σ^2 I). 

In practice, it is uncommon to know A precisely. Typical causes of uncertainty are sampling error, measurement error, human error, modeling error, or rounding error. There were attempts at addressing this model limitation in <cit.> but these relied on small magnitude errors to use a truncated series approximation for an inverse. Total least square (TLS) was developed in response to this lack of symmetry in uncertainty. Rather than isolating noise to the observed signal or right hand side, the TLS model allows for an uncertain operator and is given by (A + ) x̅ = b + η where A and b are observed with  and η a random matrix and vector, respectively. We assume that  and η have been scaled to have the same variance; see <cit.> for a modern overview of the topic. The TLS solutions, x̂_TLS, solves 

    min_x,,η   [,   η] _F   subject to    ( + )x = b + η

where [,  η] ∈^m × (n+1), A and b are the observed matrix/signal pair, and ·_F is the Frobenius norm. It was shown in <cit.> that (<ref>) can be solved via a singular value decomposition (SVD) with the closed form solution of

    x̂_TLS = (^T  - σ_n+1^2 I)^-1^T b

where σ_n+1∈ being the smallest singular value of the augmented matrix [A,  b] ∈^m × (n+1). Similar to OLS, the TLS solution yields a MLE for the model of Gaussian noise in A and b. It should be noted that (^T  - σ_n+1^2 I) is necessarily worse conditioned than ^T. Since ^T is positive definite by virtue of  being full rank, all eigenvalues of ^T are shifted closer to zero by the amount of σ_n+1^2. For a small spectral gap between σ_n and σ_n+1, the matrix will be close to singular making solutions extremely sensitive to perturbations in A. 

This can be understood intuitively; uncertainty in A permits additional degrees of freedom allowing the model to “fit” noise. To illustrate, consider the simple linear regression problem. We have access to several regressor / response pairs (a, b) ∈^2. Suppose our data is generated from the model b = a · 0 + δ with δ∼ Uniform({-1,1}) (the zero function plus discrete and uniform noise). We'd like to recover the true slope parameter, x= 0. In this instance, our three samples are given by (-0.10, 1.00), (0.00, -1.00), and (0.11, 1.00). The OLS solution returns a slope of x̂_OLS≈ 0.45 whereas TLS gives x̂_TLS≈ 297.79. This is shown graphically in Figure <ref>. TLS's ability to vary the operator results in extreme sensitivity to data provided. Although our example appears exceptional, this behavior is typical of TLS. Because of this sensitivity, a number of alternatives have been studied for uncertain operators as mentioned above. A more traditional approach for addressing the ill-conditioning encountered in TLS is ridge regression or Tikhonov regularization which we consider below.





 §.§ The Trouble with Tikhonov Regularization
 
Poor conditioning as observed in TLS is often combated with Tikhonov regularization or ridge regression (RR) whereby solutions with large norms are penalized via an ℓ_2 regularization term. The solution to the ridge regression problem solves

    min_x { x-b^2 + λ x^2 }

with λ≥ 0. The minimizing x has a closed form solution of 

    x̂_RR_λ = (^T  + λ^2 I)^-1^T b.

There is a structural connection between the RR and TLS; their solutions are nearly identical with TLS subtracting from the diagonal of ^T and RR adding to it. This can be understood from a Bayesian statistics point of view. Large values of λ correspond to stronger evidence of lower variance in x (and zero mean). In the case of an uncertain data matrix, we should have less confidence in low variance. 

Golub, Hansen, and O'Leary explored the interplay between TLS and RR in <cit.>. They considered a generalized version of the regularized total least squares (RRTLS) problem 

    min_x,,η      [,   η] _F  
    subject to    ( + )x = b + η  and  x≤γ.

The second constraint is equivalent to imposing a regularization term on the objective of the TLS problem. When the inequality constraint is replaced by an equality, the solution to (<ref>), x̂_RRTLS_α say, is given by the x that solves

    x = (A^T A + α I  )^-1 A^T b   with  α = μ (1+γ^2) - b-A  x  ^2/1+γ^2

and μ the corresponding Lagrange multiplier. Note that for all γ >0, we have α≥ -σ_n+1^2. Similarly, α = -σ_n+1^2 when γ = x̂_TLS and α = 0 when γ = x̂_OLS. This suggests that regularized total least squares merely adds a standard regularization term to the poorly conditioned TLS matrix, A^T A - σ_n+1^2 I.

Given the interplay between TLS, RR, and RRTLS, it is reasonable to solve either TLS or RR since they are effectively the same problem with different parameters; the additional term becomes regularizing when α = λ^2 > 0. Although TLS is most appropriate for the case of an uncertain matrix, its “deregularizing” effect inflames issues with conditioning making it an unpopular choice for solving typical linear inverse problems. 

There is also the tricky consideration of choosing an appropriate regularization parameter, whether in the standard RR form of (<ref>) or γ in (<ref>). Although there are many approaches for choosing a parameter such as Morozov's discrepancy principle (MDP), the unbiased predictive risk estimator method (UPR), the generalized cross validation method (GCV), or the “elbow” or “L” method to name a few, parameter selection is based on heuristic arguments or requires unknown information a priori. A detailed treatment of the above methods and their analysis can be found in <cit.>. Through the remainder of this paper, we drop our discussion of RRTLS, focusing instead on TLS and RR.




 §.§ Robust Least Squares

The central focus of robust optimization (RO) is to find a solution that is feasible over all possible realizations of uncertain variables. In our case, ∈ is unknown where = {M ∈^m × n: |M| ≤ D } and D is an element-wise positive matrix chosen in a principled fashion, i.e., the degree to which matrix entries are quantized. Our robustified version is written as a minimax optimization problem in (<ref>) with the appropriate . The inner maximization problem guards against over-fitting, effectively regularizing our solution, and reducing sensitivity to perturbations in .
Note that RO avoids placing a statistical prior on Δ, which can be a strength or weakness depending on the model.

The outline of the rest of the paper is as follows: in section II we derive a closed form solution to the inner maximization problem thereby showing its tractability, section III discusses computational methods for solving problem (<ref>), and section IV provides the results to numerical experiments.



§ CLOSED FORM SOLUTION OF INNER OBJECTIVE AND THEORY



 §.§ Floating Point Uncertainty

We begin by breaking (<ref>) into an inner maximization and outer minimization problem:

    min_x  {max_ || ≤ D   (+)x -b ^2 }_f(x)

which is equivalent to

    min_x   f(x)   subject to    f(x) = max_|| ≤ D(A+)x - b ^2.

Here || ≤ D indicates that the magnitude of elements in  are bound by the corresponding non-negative components of D. That is, D_i,j constrains element _i,j of the uncertainty matrix. The inner maximization problem will recover f(x). Because the function x ↦ (+)x -b ^2 is convex in x and the supremum over an arbitrary family of convex functions is convex, it follows that f is convex making min_x f(x) an unconstrained convex optimization problem. To evaluate f and find a subgradient, it remains to find the maximizing ∈^m × n. 
 
	The maximizing function f(x) in (<ref>) is given by 
	
    f(x) = Ax-b ^2 + 2 ⟨  |Ax-b|, D|x|  ⟩ +   D|x|  ^2

	where |x| and |Ax-b| denote the vectors of component-wise absolute values. 
	
		Since f(x) = max_|| ≤ D(A+)x - b ^2 must be maximized over Δ, we can fix x, define c:=  x - b, then treat it as constant. Exploiting row-wise separability, we write 
		
    max_|| ≤ D  x + c ^2 = ∑_i = 1^m max_|^T| ≤ D_i^T  (^T x + c_i)^2
  
		where ^T  and D_i^T are the i row of  and D, respectively, and c_i is the i element of vector c. We now work row by row. Note that we can switch to absolute values rather than squares when maximizing for each row. Applying the triangle inequality gives an upper bound
		
    |^T x + c_i|  ≤      |c_i| +  ∑_j = 1^n |_i,j|   | x_j| 
    ≤     |c_i| + ∑_j = 1^n D_i,j |x_j|    =    D_i^T 
    		|x|  + |c_i|.

		It is easily verified that the upper bound is achieved when
		
    = D ⊙(x   c^T)

		making it a solution to the inner maximization problem. Here, ⊙ denotes the Hadamard or element-wise product. Recalling that c_i = [Ax - b]_i, we simplify to 
		
    f(x) =  Ax-b ^2 + 2 ⟨  |Ax-b|, D|x|  ⟩ +   D|x|  ^2.

	

Using Theorem <ref>, the optimization problem in (<ref>) can be rewritten as 		

    min_x f(x) = min_x {Ax-b ^2 + 2 ⟨  |Ax-b|, D|x|  ⟩ +   D|x|  ^2      }.

In general, f is not differentiable (e.g., let m=n=1 and A=D=b=1, then f(x)=(x-1)^2 + 2|x|·|x-1| + x^2 is not differentiable at x∈{0,1}), but we are guaranteed a subgradient by virtue of convexity. Furthermore, a generalization of Danskin's Theorem <cit.> provides us with a method to find elements of the subgradient for all x. In particular,

    f'(x) := 2 (+_x)^T[(+_x)x - b ]  ∈ ∂ f(x)

where _x indicates the optimal  for a given x as provided in Eq. (<ref>).



 §.§ Fixed Point Uncertainty

Fixed point uncertainty is a special case of (<ref>) with all elements bound by the same value. We can write this constraint as Δ_∞≤δ with _∞ representing the largest magnitude element of . Problem (<ref>) becomes

    min_x  {max_  (+)x -b ^2 }

with a solution denoted by x̂_RO. We focus on fixed point error for the remainder of the paper. This instance gives rise to the following corollary.


	For fixed point uncertainty, problem (<ref>) reduces to 
	
    min_x { x-b^2 + 2 δx_1  x-b_1 + m δ^2 x_1^2 }.

	
		We take _k ∈^k to be the vector composed of ones. Note that Δ_∞≤δ is equivalent to |Δ| ≤ D when D =  δ _m^_n^T. By Theorem <ref>, the result follows easily using properties of inner products.
	




 §.§ Explicit Regularization of Robust Objective

The robust formulation presented above is intended to address uncertainty in the data matrix from rounding error. Although implicit regularization occurs via the 3rd term in both (<ref>) and (<ref>), we don't address poor conditioning of A directly. To illustrate this point, note that x̂_RO→x̂_OLS as δ→ 0. If A is poorly conditioned, such behavior is undesirable. This can be remedied by adding an ℓ_2 regularization term and solving


    min_x {(  x-b^2 + 2 δx_1  x-b_1 + m δ^2 x_1^2 ) + λ^2 x^2      }.


Since the robust objective and regularization term are both convex, so is their sum. Furthermore, the modified objective is 2 λ^2 strongly convex yielding a unique minimizer, x̂_RRO_λ. Strong convexity also implies a bounded sequence of iterates, i.e., x_k - x^*≤ M,  and therefore guarantees convergence under the mirror-descent method. In fact, the same technique works with regularizers other than ℓ_2, and can be solved via proximal mirror-descent methods when the regularizer is smooth or has an easy-to-compute proximity operator <cit.>.



§ ALGORITHMS
 
Equipped with an element of the subdifferential, we can employ a variety of solvers. Bundle methods are a promising choice as they sequentially form an approximation to the objective. This is done by using the location, function value, and subgradient of previous iterates to lower bound the objective with supporting hyperplanes. At each step, a direction finding quadratic program must be solved, but can be dealt with rapidly thanks to software such as CVXGEN <cit.> and ECOS <cit.>. A survey on the topic can be found in <cit.>.

Since the problem is unconstrained, another option is to use smooth optimization techniques, and in particular quasi-Newton methods that form low-rank approximations of ∇^2 f. Our objective f is not differentiable, much less twice so, and consequently methods requiring second derivatives seem theoretically inappropriate. However, there is a growing body of literature going back to Lemaréchal <cit.> recognizing the empirical success of these methods. See <cit.> and references therein. The success of these methods depends on the smoothness near the solution. Using the  MATLAB package <cit.> with random and zero vector initializations, and using the package's default solver of limited-memory BFGS <cit.>, we observed fast and accurate convergence for several test problems. We used this method for the numerical experiments presented in the next section.

Subgradient descent is an appealing option due to its simplicity and flexibility. The method can easily handle constraints via projections or a regularization term h(x) with proximal operators. Recent advances in proximal subgradient descent methods and their analysis can be found in <cit.> and <cit.>. Convergence may be slow, but is less of a concern with heavy quantization. The 3rd terms in (<ref>) and (<ref>) effectively regularize the objective when δ and D are large yielding faster convergence. A drawback of the subgradient descent method is that a step-size scheme must be chosen in advance. One popular choice is the diminishing, non-summable step length given by t_k = 1/√(k+1)·f'(x_k). A discussion on the choice of step lengths can be found in <cit.>. Stopping criteria are also difficult to specify since it is not a true descent method (the objective is not guaranteed to decrease on every iteration). The algorithm therefore typically runs for a fixed number of iterations.



§ NUMERICAL EXPERIMENTS
 


  
*SetupWe use MATLAB's pseudo-random number generator  to draw 10,000 standard normal matrices, A̅ ∈ ^30 × 15. The bar notation in this section indicates true and unobserved values. We fixed the condition number of our matrices at 100. We did so by performing a singular value decomposition (SVD) such that A̅ = U Σ V^T where U, V are unitary and Σ is diagonal, then replaced the diagonal of Σ by linearly decaying values from 1 →1/100. This ensures our test matrices show sensitivity to noisy measurements prior to quantization.



We conducted two experiments. In the first, we draw random vectors from a heavy-tailed Cauchy distribution with median 0 and scale parameter 1. In the second, the true solution x̅ has its first element drawn from x̅_1 = ± 100 with a sign drawn uniformly at random and the remaining elements drawn from a standard normal distribution, x̅_2:15∼𝒩(0,I). This setup exposes the bias of RR solutions. We remark that we do this for illustrative purposes; if one suspected a signal of having such a large element, it would of course make sense to run an outlier detection method first.


We obtain our “true” right hand side by taking the image of x̅ under A̅, that is, b̅ := A̅x̅, then generate our observed measurement vector by letting b := b̅ + η where η∼𝒩(0, b̅^2/m ·SNR I) with the desired signal-to-noise ratio (SNR) fixed at SNR = 50. Finally, we quantize A̅ which yields our observed A, then solve problems for OLS (<ref>), TLS (<ref>), RR_λ (<ref>), RO (<ref>), and RRO_λ (<ref>). The λ for RR_λ and RRO_λ are chosen according to the GCV and MDP criteria discussed in Section <ref>, and the δ for RO and RRO_λ is based on the observed accuracy, e.g., if A̅ is rounded to the 2^nd decimal, then  δ = 0.5 · 10^-2. We use the  implementation of limited-memory MFGS to solve the robust problems with a random initialization. The GCV parameters is recovered using an approximate line search for the corresponding objectives provided in Vogel's text <cit.>. We omit UPR because of its similarity to the GCV parameter.

The MDP parameter is found using the bisection method to solve A  x̂_RR_λ - b^2 - ρ = 0 for λ, where x̂_RR_λ is the solution to (<ref>) and the parameter ρ ideally reflects the residual's noise floor. We chose this method since it is computationally cheap and easy to program, though it is also possible to solve a constrained least-squares problem directly using standard software like  <cit.>. 
The goal is to find λ such that A x̂_RR_λ- b^2 equals the expected uncertainty introduced by noise.
For A̅x̅ - b = η∼𝒩(0, σ^2 I) which is an m component Gaussian random variable, 

    1/σ^2η^2 = ∑_i=1^m ( η_i/σ)^2 ∼χ^2(m)
 
since η_i/σ∼𝒩(0,1) and ( η_i/σ)^2 ∼χ^2(1). When m=30, as in our case, both the mean and median are approximately 30 (exact for mean). To ensure feasibility, i.e., a real root, we choose ρ such that ℙ(A̅x̅-b^2 - ρ≤ 0 ) = 95%, which can be easily calculated using χ^2 inverse CDF tables. For our experiments, ρ≈2 b̅^2/3 ·SNR. An appropriate ρ is generally unknown a priori unless one knows the noise variance.



  
*Numerical Results
Results for the Cauchy distributed measurement vector are displayed in Figure <ref>, which shows mean error as a function of rounding digit for A̅, and Figure <ref>, which illustrates an empirical probability density function (pdf) of component-wise errors for different methods using a standard Gaussian kernel. In this case, A is quantized to the hundredth spot. Note that the robust solution does well at reducing error in the face of heavy quantization. This is evidenced by low mean error through the thousandth spot and small variance as seen in the pdf.
The RO solution converges to the OLS solution as expected when δ→ 0, i.e., when the quantization effect is small. This might be undesirable, especially when A is ill-conditioned, and can be addressed via the inclusion of a regularization term. The mean relative error performance for the regularized robust problem given by (<ref>) is shown in Figure <ref>. Since the proposed method does not improve estimation universally, we recommend its use when quantization results in precision loss of greater than 0.1%. We reiterate that when quantization effects are heavy, there is a notable benefit of the proposed method and the parameter δ can be chosen in a natural way. 




A second observation is that RO doesn't sacrifice accuracy for bias as RR does. This phenomenon can be observed in Figure <ref> which depicts empirical pdf's of component- wise error for different methods. The plots to the right show error behavior for small components of the solution, that is, for x̅_2:15∼𝒩(0,I). Since values are close to zero, RR is expected to perform well at estimating the true solution. Indeed, RR implicitly assumes a prior with mean zero. 
The densities on the left show errors for “large” components drawn from {± 100} (and the sign of the error is adjusted by the sign of the large component, so a negative error indicates that the estimate is biased toward zero, and a positive error indicates bias to ±∞). Rather than localizing about zero, the absolute error's mode is observed around -7 for GCV and UPR and -25 for MDP. Heavier penalty terms (bigger λ) place less weight on the LS term in the objective and bias estimates more aggressively towards zero. Typical λ values for GCV, UPR, and MDP are 0.041, 0.042, and 0.17, respectively. The robust solution accurately estimates the large and small components values without the extreme errors observed with OLS and TLS; that is, RO appears to have a low-variance, in-line with RR, but also low-bias, comparable or better than OLS.




§ CONCLUSION

In this paper, we presented a robust method for addressing fixed and floating point uncertainty due to quantization error. After reviewing the limitations of existing methods, we formulated a tractable robust objective and presented algorithms to solve it. The only parameter necessary in our formulation is chosen in a principled fashion, i.e., by observing the degree to which matrix elements are rounded. Our numerical experiments show that robust least squares outperforms OLS, TLS, and RR, effectively balancing bias and accuracy. 



 §.§ Acknowledgments

We would like to thank Marek Petrik for his keen eye and helping us spot an elegant simplification to our approach.  




ieeetr




