SLIDE 1
Condition estimation of linear algebraic equations and its application to feature selection
Joab Winkler Department of Computer Science, The University of Sheffield, Sheffield, United Kingdom
The Institute of High Performance Computing A* Agency for Science, Technology and Research Singapore January 2019
SLIDE 2 1
Introduction Mathematical background
2
Regression
3
Condition numbers and regularisation The effective condition number The discrete Picard condition Tikhonov regularisation Componentwise condition numbers
4
Feature selection Feature selection and condition estimation
5
Summary
SLIDE 3
Introduction
Several problems require the prediction of the output of a physical system for which the sample size n is much smaller than the dimension of the data p: Chemometrics Brain imaging Genomics Gene selection from microarray data Text analysis The condition n < p implies that there are many models that satisfy the given data and important issues therefore arise: Which model from this infinite set of models should be chosen? What is the criterion that should be used for this selection? Can the selection be generic, that is, not problem dependent, such that prior information is not required?
SLIDE 4 Mathematical background
These problems yield an equation of the form Ax = b + ε, A ∈ Rm×n, b ∈ Rm, x ∈ Rn where m < n, rank A = m and ε is the noise. The least squares minimisation of ε leads to the normal equation ATAx = ATb whose solution is xsoln = A†b = VS†UTb where the superscript † denotes the pseudo-inverse. The solution is xsoln = xln + x0, xln = V
1 UTb
0n−m
x0 = V 0m r
- where xln is the minimum norm solution, x0 lies in the null space of A, r
is arbitrary, and S = S1 0m,n−m
SLIDE 5 The solution xln is unsatisfactory for two reasons: Prediction accuracy: This solution may have low bias and high
- variance. Prediction accuracy can sometimes be improved by
reducing or setting to zero some coefficients of xln. Interpretation: It is usually desirable to choose the most important components of xln that characterise the physical system being considered. Methods that are used to overcome these problems: Ridge regression: The magnitude of the components of xln is reduced continuously:
It is more stable than subset selection. It does not set any components to zero and thus it does not yield a sparse model that can be easily interpreted.
Subset selection: Components of xln are deleted in discrete steps:
The models are strongly dependent on the components that are deleted because the elimination procedure is discrete. A small change in the data can cause a large change in the selected model, which reduces the prediction accuracy.
SLIDE 6 Ridge regression (Tikhonov regularisation) The sensitivity of the solution xln to perturbations in b can be reduced by a constraint on the magnitude of the solution xreg: xreg = arg min
x
- (Ax − b)T (Ax − b) + λ x2
, λ > 0 and thus
λ > 0 The lasso (‘least absolute shrinkage and selection operator’) The retension of the advantages of ridge regression (stability) and subset selection (sparsity) are combined in the lasso: xlasso = arg min
x
- (Ax − b)T (Ax − b)
- subject to
x1 ≤ t which can also be written as xlasso = arg min
x
- (Ax − b)T (Ax − b) + λ x1
- ,
λ > 0
SLIDE 7 The elastic net This method is an improvement on the lasso and it combines L1 and L2 regularisation: xelastic = arg min
x
- (Ax − b)T (Ax − b) + λ1 x1 + λ2 x2
where λ1, λ2 > 0 The solutions from Tikhonov regularisation, the lasso and the elastic net reduce the sensitivity of the least norm solution xln to perturbations in b, but there are differences between these forms of regularisation.
SLIDE 8 Compare Tikhonov regularisation Tikhonov regularisation imposes a Gaussian prior on the parameters
Tikhonov regularisation does not impose sparsity on xreg. The solution xreg has a closed form expression. with the lasso The lasso imposes a Laplacian prior on the parameters of the model. The lasso favours sparse solutions because some coefficients of xlasso are set to zero. The sparsity of xlasso increases as λ increases. The solution xlasso does not have a closed form expression and quadratic programming is required for its computation. and the elastic net The sparsity of xelastic is similar to the sparsity of xlasso. The solution xelastic favours a model in which strongly correlated predictors are usually either all included, or all excluded. The solution xelastic is much better than xlasso for some problems.
SLIDE 9
A regularised solution (Tikhonov, the lasso and the elastic net) is stable with respect to perturbations in b, but several points arise: Is regularisation always required when the data b are corrupted by noise? Must specific conditions on A and b be satisfied in order that regularisation is imposed only when it is required? What are consequences of applying regularisation when it is not required? If regularisation is required, then rmethod = xln − xmethod = 0, method = {reg, lasso, elastic} Can bounds be imposed on rreg, rlasso and relastic, such that these errors induced by regularisation are quantified? The answers to these questions are most easily obtained if Tikhonov regularisation is considered because the constraint in the 2-norm lends itself naturally to the SVD.
SLIDE 10
Regression
The use of regularisation is usually justified for three reasons: It reduces or eliminates over-fitting in regression. It reduces the sensitivity of the regression curve to noise in the data. It imposes a unique solution in feature selection Ax = b where A ∈ Rm×n, b ∈ Rm, x ∈ Rn, m < n and rank A = m. But There are well-defined problems for which regularisation must not be used because it causes a large degradation in the solution. and thus Can a quantitative test be established such that regularisation is used only when it is required? Regression provides a good example of the correct use, and the incorrect use, of regularisation.
SLIDE 11 Example 1 Consider the points (xi, yi) , i = 1, . . . , 100, where the independent variables xi are not uniformly distributed in the interval I = [1, . . . , 20], the dependent variables yi are given by yi =
33
ak exp
2σ2
d
i = 1, . . . , 100 the centres dk of the 33 basis functions are uniformly distributed in I and σd = 1.35. Consider two sets of data points, y = y1 and y = y2, and the perturbations δy1 and δy2, δy1, δy2 ∼ N
and δy1 y1 = 3.41 × 10−6 and δy2 y2 = 8.27 × 10−4
SLIDE 12 10 20 30 40 50 60 70 80 90 100 i
50 100 150 200 250
5 10 15 20 25 30 35 i
0.5 1 104 exact data perturbed data
Figure: The exact curve (top), and the coefficients ai (bottom) for the exact data y = y1 and the perturbed data y = y1 + δy1.
SLIDE 13 10 20 30 40 50 60 70 80 90 100 i
0.5 1 5 10 15 20 25 30 i
0.5 1 107
Figure: The exact curve (top), and the coefficients log10 |ai| (bottom) for the exact data y = y2 and the perturbed data y = y2 + δy2.
SLIDE 14
The interpolated curve is unstable for the data set y = y1. The interpolated curve is stable for the data set y = y2. The coefficient matrix A ∈ R100×33 is the same for y = y1 and y = y2, and its condition number is κ(A) = 5.11 × 108 Thus The presence of noise in the vector b, where Ax = b, does not imply that x is sensitive to changes in b. The condition κ(A) ≫ 1 does not imply that the equation Ax = b is ill-conditioned. Tikhonov regularisation yields a very good result (numerically stable and a small error between the theoretically exact solution and the regularised solution) for y = y1, but an unsatisfactory result for y = y2 (a very large error between the theoretically exact solution and the regularised solution).
SLIDE 15
Condition numbers and regularisation
The 2-norm condition number of A ∈ Rm×n is κ(A) = s1 sp , p = min(m, n) where si, i = 1, . . . , p, are the singular values of A and rank A = p. The condition number κ(A) cannot be a measure of the stability of Ax = b because it is independent of b. It is necessary to develop a measure of stability that is a function of A and b. This leads to:
A refined normwise condition number - the effective condition number - which is a function of A and b. Componentwise condition numbers - one condition number for each component of x.
SLIDE 16 The effective condition number
The effective condition number η(A, b) of ATAx = ATb, A ∈ Rm×n, m ≥ n is a refined normwise condition number. Theorem 1 Let the relative errors ∆x and ∆b be ∆x = δx x and ∆b = δb b The effective condition number η(A, b) of AT Ax = ATb is equal to the maximum value of the ratio of ∆x to ∆b with respect to all perturbations δb ∈ Rm η(A, b) = max
δb∈Rm
∆x ∆b = 1 sn c S†c where A = USV T is the SVD of A and c = UTb.
SLIDE 17 Proof It follows from ATAx = ATb that x = VS†UTb = VS†c and thus δx =
sn and the division of both sides of this inequality by x yields δx x ≤ 1 sn δc c c x = 1 sn δb b c x and thus η(A, b) = max
δb∈Rm
∆x ∆b = 1 sn c S†c = max
δb∈Rm
∆x ∆b = 1 sn b S†UTb
SLIDE 18 η(A, b) = max
δb∈Rm
∆x ∆b = 1 sn c S†c = max
δb∈Rm
∆x ∆b = 1 sn b S†UTb The minimum value of η(A, b) occurs when c = en (ei is the ith unit basis vector). If some conditions on b are satisfied, the maximum value of η(A, b)
- ccurs when c = e1, and thus
1 ≤ η(A, b) ≤ s1 sn = κ(A) More generally, since b = Uc: η(A, b) ≈ 1 if the dominant components of b lie along the columns ui of U that are defined by large values of i. η(A, b) ≈ κ(A) if the dominant components of b lie along the columns ui of U that are defined by small values of i.
SLIDE 19 The discrete Picard condition
η(A, b) = max
δb∈Rm
∆x ∆b = 1 sn c S†c, c = UTb Consider the ratio |ci |/si, i = 1, . . . , n.
i
(c) (b) (a)
|ci|/si
Figure: Three forms for the ratio |ci |/si: (a) monotonically decreasing, (b) constant and (c) monotonically increasing.
SLIDE 20
If |ci|/si → 0 as i → n, then the constants |ci| decay to zero faster than the singular values si. The condition |ci| si → 0 as i → n is called the discrete Picard condition. If this condition is satisfied, then η(A, b) ≈ κ(A) and thus Ax = b is ill-conditioned. Tikhonov regularisation can be used to yield a well-conditioned solution x. If |ci|/si ≈ 1, i = 1, . . . , n, then η(A, b) ≈ κ(A)/√n. Tikhonov regularisation cannot be used to yield a well-conditioned solution x. If |ci+1| ≫ |ci| , i = 1, . . . , n − 1 then η(A, b) ≈ 1 and regularisation must not be applied.
SLIDE 21 Computation of the discrete Picard condition The important term for the evaluation of the discrete Picard condition is x =
- A†b
- =
- S†c
- If noise is present, then the square of the magnitude of the perturbed
solution is x + δx2 =
- A†(b + δb)
- 2 =
- S†(c + δc)
- 2 =
n
ci + δci si 2 Assume the magnitude of the perturbation |δci| is approximately constant, such that |δci| ≪ |ci| , i = 1, . . . , p − 1 |δci| ≈ |cp| , i = p |δci| ≫ |ci| , i = p + 1, . . . , n
SLIDE 22 It follows that if the discrete Picard condition is satisfied |ci + δci| si ≈
|ci | si ,
i = 1, . . . , p − 1
|cp+sp| sp
, i = p
|δci | si ,
i = p + 1, . . . , n
i
(c) (b) (a)
Figure: The ratios (a) |ci |/si, (b) |δci |/si and (c) |ci +δci |/si if the discrete Picard condition is satisfied.
SLIDE 23
If the discrete Picard condition |ci| si → 0 as i → n is satisfied, the function |ci+δci |/si cannot be computed because it is sensitive to noise. A similar result occurs if |ci| ≈ si, i = 1, . . . , n because |ci + δci| si ≈ |δci| si , i = p + 1, . . . , n If |ci+1| ≫ |ci| , i = 1, . . . , n − 1 computational problems do not occur because |ci + δci| si ≈ |ci| si , i = 1, . . . , n
SLIDE 24
Summary (discrete Picard condition) The form of the term |ci|/si defines the stability of Ax = b. If |ci| si → 0 as i → n then Ax = b is ill-conditioned. The dominant components of b lie along the columns ui of U that are defined by small values of i (large singular values). Tikhonov regularisation enables a well-conditioned solution to be computed. If |ci+1| ≫ |ci| , i = 1, . . . , n − 1 then Ax = b is well-conditioned and |ci|/si is numerically stable. The dominant components of b lie along the columns ui of U that are defined by large values of i (small singular values).
SLIDE 25 Example 2 Consider the regression problem in Example 1. The data points y = y1 yield an ill-conditioned equation AT Ax = ATb. The effective condition number is η(A, b) = 4.61 × 108. The data points y = y2 yield a well-conditioned equation AT Ax = ATb. The effective condition number is η(A, b) = 7.94.
5 10 15 20 25 30 i
perturbed data exact data
Figure: The ratio log10 |ci |/si for the exact data y = y1 and the ratio log10 |ci +δci |/si for the perturbed data y = y1 + δy1.
5 10 15 20 25 30 i 1 2 3 4 5 6 7
Figure: The ratio log10 |ci |/si for the exact data y = y2 and the ratio log10 |ci +δci |/si for the perturbed data y = y2 + δy2.
SLIDE 26 Tikhonov regularisation
If the discrete Picard condition is satisfied, the equation AT Ax = ATb is ill-conditioned and Tikhonov regularisation requires the solution x(λ) of (AT A + λI)x(λ) = ATb and thus x(λ) = V
UTb, λ ≥ 0 What is the error between the regularised solution (λ > 0) and the exact solution (λ = 0) if:
the discrete Picard condition is satisfied the discrete Picard condition is not satisfied
Assume that λ∗ is the optimal value of the regularisation parameter λ.
SLIDE 27 If the discrete Picard condition, |ci|
si → 0 as i → n, is satisfied, the
solution x(0) is dominated by the large singular values and it is independent of the small singular values. The error is small x(λ∗) − x(0) x(0) ≈ λ∗ λ∗ + s2
1
≪ 1 because Tikhonov regularisation filters out the small singular values. If |ci| ≈ si, i = 1, . . . , n, the error is smaller because all the singular values contribute to x(0) x(λ∗) − x(0) x(0) ≈ n − p n 1
2
Tikhonov regularisation filters out the components of x(0) associated with the small singular values, but the components associated with the large singular values are not affected. If |ci+1|/si+1 ≫ |ci|/si, i = 1, . . . , n − 1, the error is large x(λ∗) − x(0) x(0) ≈ 1, x(λ∗) ≈ 0 because x(0) is dominated by the small singular values, all of which are filtered out by Tikhonov regularisation.
SLIDE 28 Componentwise condition numbers
The condition number κ(A) and effective condition number η(A, b) are defined in the normwise sense. Condition numbers defined in the componentwise sense are more refined. A condition number is assigned to each component of x, where AT Ax = ATb. The condition number ρ(A, b, xj) of the component xj of x is defined as ρ(A, b, xj) = max
δb∈Rm
∆xj ∆b , j = 1, . . . , n where ∆xj = |δxj| |xj| and ∆b = δb b
SLIDE 29 Theorem 2 The condition number ρ(A, b, xj) of the component xj, j = 1, . . . , n, of AT Ax = ATb is ρ(A, b, xj) =
j A†
b |xj| =
j A†
b
j A†b
1 cos γj > 1 where γj is the angle between b and the jth row of A†. Proof If ej is the jth unit basis vector, then a change δb in b causes a change δxj in xj δxj = eT
j δx = eT j A†δb,
j = 1, . . . , n and thus ∆xj = |δxj| |xj| ≤
j A†
δb |xj| , j = 1, . . . , n and the result follows.
SLIDE 30 Example 3 Consider the regression problem in Example 1. The data points y = y1 yield an ill-conditioned equation AT Ax = ATb. The data points y = y2 yield a well-conditioned equation AT Ax = ATb.
5 10 15 20 25 30 35 i 6 7 8 9 10 11 effective condition number componentwise condition numbers
Figure: The effective condition number and the componentwise condition numbers, on a logarithmic scale, for y = y1.
5 10 15 20 25 30 35 i 5 10 15 20 componentwise condition numbers effective condition number
Figure: The effective condition number and the componentwise condition numbers for y = y2.
SLIDE 31 Feature selection
Many problems in feature selection yield the equation Ax = b + ε where A ∈ Rm×n, b ∈ Rm, x ∈ Rn, m < n, rank A = m and ε is the noise in the data b. The normal equations are
This equation has an infinite number of solutions xsoln xsoln = V
1 UTb
0n−m
r
where xln is the solution of minimum norm, x0 lies in the null space of A and is orthogonal to xln, r ∈ Rn−m is arbitrary and A = USV T = U
0m,n−m
S1 = diag {si}m
i=1
SLIDE 32 The most important features (components of xsoln) of the system are usually sought and a sparse solution is therefore desired. Use the n − m degrees of freedom in r in the null space vector x0 to
- btain a sparse and regularised solution xsoln
xsoln = V
1 UTb
0n−m
r
The solutions xlasso and xelastic ignore the vector x0. The vector r is chosen such that if the d components k1, k2, . . . , kd,
- f xln have large condition numbers, then these d components of the
solution xsoln are equal to zero. The lasso and the elastic net yield sparse solutions from constrained minimisation problems but they are approximate solutions of the normal equations
The need, or otherwise, for regularisation, and the errors between (a) the solutions xsoln and xlasso, and (b) the solutions xsoln and xelastic, are not considered.
SLIDE 33 Feature selection and condition estimation
Since xsoln = xln + x0 = xln + V
r
V1 V2 0m r
xsoln = xln + V2r, V2 ∈ Rn×(n−m) By definition, the components of xln that have the d largest condition numbers, 1 ≤ d ≤ n − m, satisfy xln,k + x0,k = 0, k = k1, k2, . . . , kd that is, these components of xsoln are set to zero, thereby inducing sparsity in xsoln. The vector r is chosen to satisfy these d equations.
SLIDE 34 Summary
A large condition number of A, κ(A) ≫ 1, does not imply that the equation Ax = b is ill-conditioned. The effective condition number η(A, b) is a better measure of the stability of Ax = b because it is a function of A and b, which must be compared with κ(A), which is a function of A only. If the discrete Picard condition is satisfied, then η(A, b) ≈ κ(A) and Tikhonov regularisation yields a stable solution whose error is small. More refined measures of the stability of Ax = b are the componentwise condition numbers, which measure the stability of each component of x due to a perturbation in b. The equation Ax = b for feature selection has an infinite number of
- solutions. The lasso and elastic net apply regularisation to the least
norm solution xln, and they do not consider the vectors xnull that lie in the null space of A. Consideration of the componentwise condition numbers of xsoln = xln + xnull yields a sparse solution that satisfies the normal equations, and regularisation is satisfied in the componentwise sense.
