"There are three kinds of lies: lies, damned lies, and statistics." The origin of the statement is lost in time, but the sentiment lives on. I'm tempted to give it a 21st century update. "There are three kinds of lies: lies, damned lies, and floating point arithmetic."
I am currently working on a research project in which the following problem arises. We start with an $m \times n$ matrix $A$ with $m < n$, and a set of indices $I\subseteq \{1,\dots,n\}$. I need to confirm that $x_i = 0$ for all $i\in I$ and for all $x\in \mathcal{N}$, where $\mathcal{N} = \{x\in \mathbb{R}^n : Ax = 0\}$. There are various names for $\mathcal{N}$; I'll go with the null space of $A$.
How do you test whether a particular component of every vector in an uncountably infinite set is zero? Since $\mathcal{N}$ is a vector space, we can compute a basis for it and just check the basis vectors. If $x_i =0$ for every basis vector $x$ of $\mathcal{N}$, it will be true for every vector in $\mathcal{N}$.
This is going on in a Java program, and I have tried multiple linear algebra packages. Initially I had the idea to use a QR decomposition of $A$ to get a basis for the null space, but that requires (at least as far as I know) an upper triangular $R$ component, and neither of the two most promising packages I tried delivered one. (I'm pretty sure the Javadocs for one or both said $R$ would be upper triangular ... but it wasn't.) So I switched to singular value decomposition (SVD) using the EJML library. SVD is generally slower than QR decomposition, but tends to be more accurate.
I got an answer, but a check performed by the Java program warned me that some of the $x_i$ values that were supposed to be zero were in the vicinity of $10^{-6}$, which is a bit dicey to shrug off as rounding error. (The basis vectors are all normalized, by the way.) So I decided to port $A$ and some other stuff over to R and do some testing there ... which is when the adventure began.
The dimensions of this $A$ matrix were $m=401$ and $n=411$. EJML gave me a basis of 17 vectors, indicating that $A$ had nullity 17 and thus rank 394. (The nullity of matrix $A$ is the dimension of $\mathcal{N}$, and is in turn the difference between $n$ and the rank of $A$.) I tested this by doing QR and SVD decompositions of $A$ (using the qr() and svd() functions from the R base package) and the Rank() function from the R pracma package. The QR and SVD results indicated that $A$ had rank 393 and nullity 18. Only pracma::Rank() agreed with EJML ... and that came with the warning "Rank calculation may be problematic". (As if the rest of this stuff isn't.)
So EJML might have shorted me one dimension in the null space ... or maybe not. The smallest eigenvalues of $A$ from the SVD were as follows: 9.16e-05, 5.37e-06, 1.09e-10, 1.13e-15, 7.56e-16, 6.20e-16, 4.87e-16, 3.79e-16, 3.37e-16, 3.16e-16. So the last seven (italicized) were pretty clearly zero and the preceding one (bold) is maybe zero, maybe not. The SVD does not produce the ten zero eigenvalues resulting from the difference between $m$ and $n$, so basically we're looking at a nullity of 18 if you treat the bold one as zero and 17 if you don't.
The 17 basis vectors EJML gave me were all at least close to zero in all components $i\in I$. R had almost all the designated components $10^{-7}$ or smaller. The pracma::nullspace() function gave me 18 basis vectors (despite pracma::Rank() coming up with a nullity of 17), and some of the corresponding $x_i$ were pretty big (0.29, for example) -- way too big to be zero. So my condition is (probably) satisfied if I believe EJML and (definitely) not satisfied if I believe the R computations.
Hoping for some clarity, I calculated (in R) $Ax$ for each basis vector $x$ from either EJML or pracma::nullspace(). The products with the EJML basis vectors were all zero vectors (meaning no components larger than about $10^{-15}$ in magnitude). For the pracma basis vectors, the products had components around $10^{-9}$ in many cases. So much for clarity. Do we treat $10^{-9}$ as approximately zero and believe pracma, or do we call it nonzero and believe EJML?
I'm not new to the limitations of double-precision arithmetic. My first course in grad school was numerical analysis, and I've had a fair share of adventures with MIP models in CPLEX where the numerics were, shall we say, dicey. Still, it's rather frustrating to have to try multiple packages and then have to judge (meaning guess) which one is lying to me the least.
Showing posts with label numerical analysis. Show all posts
Showing posts with label numerical analysis. Show all posts
Tuesday, October 8, 2019
Friday, February 28, 2014
Ridge Regression Revisited
I've been neglecting the blog a bit lately, partly because I haven't had much to say and partly because I've been a bit busy with other things. One of the things keeping me occupied is the excellent Statistical Learning course offered by Stanford. (In one of those coming-full-cycle things, now that I'm retired from teaching, I'm going back to my student roots and taking MOOCs in subjects of interest.)
During one of the lectures I learned something new that was both interesting and had the side effect of making me feel old (a job my knees are entirely up to on their own, thank you very much). Ridge regression has apparently been picked up from the scrap heap of statistical history and repurposed.
To explain ridge regression as I originally learned it, let me start by defining the phrase ill-conditioned matrix. Let $A$ be a square matrix (we'll operate over the reals, although I think this works as well for complex matrices) and let $\left\Vert \cdot\right\Vert $ be a consistent matrix norm. The condition number of $A$, assuming that is nonsingular, is $$\kappa = \left\Vert A\vphantom{^{-1}}\right\Vert \left\Vert A^{-1}\right\Vert \ge 1.$$ For singular matrices, $\kappa$ is either undefined or $\infty$, according to your religious preferences.
Somewhat informally, $\kappa$ can be interpreted as a measure of how much rounding errors may be inflated then inverting $A$ (using floating-point arithmetic ... "pure" mathematicians don't commit rounding errors). An ill-conditioned or "stiff" matrix $A$ is non-singular but has a large condition number (the criterion for "large" being rather fuzzy), meaning that rounding errors during the process of inverting (or factoring, or otherwise torturing) $A$ may result in an answer with lots of incorrect digits. Long ago, a graduate school instructor told my class the story of an engineer who had inverted an ill-conditioned matrix, without realizing it was ill-conditioned, using home-brew FORTRAN code. He'd published the inverse in a journal article to four or five decimal places. All the decimal places in most of the coefficients were wrong. So were lots of digits to the left of the decimal points.
Here comes the tie-in to regression. Suppose we are trying to fit the model $$y=x'\beta+\epsilon$$(where $\epsilon$ is random noise with mean zero), and we have collected observations $(X,Y)$ with $X$ an $n\times m$ matrix ($n$ = sample size, $m$ = number of predictors). The normal equations tell us that the least squares estimate for $\beta$ is $$b=(X^\prime X)^{-1}X^\prime Y.$$In theory, $X^\prime X$ is nonsingular unless the predictors are perfectly multicollinear. In practice, if $X^\prime X$ is nonsingular but stiff, our coefficient estimate $b$ can suffer drastic rounding errors. To mitigate those errors, ridge regression tampers with $X^\prime X$, adding a small positive multiple of the identity matrix of dimension $m$. The modified normal equations are$$b_\lambda=(X^\prime X + \lambda I)^{-1}X^\prime Y.$$The result $b_\lambda$ is a deliberately biased estimate of $\beta$, in which we trade the bias for reduced rounding error. We generally choose $\lambda > 0$ as small as possible subject to the need to stabilize the numerics.
When I first encountered ridge regression, we were in the (first) mainframe era. Machine words were small, single-precision arithmetic ruled the computing universe, and rounding errors were relatively easy to generate. Somewhere along the way, computing hardware got smaller and cheaper, double precision arithmetic became the norm and separate floating point units took over the computational load. At the same time, matrix factorization replaced matrix inversion, and factoring algorithms got smarter about managing rounding error. The problems with stiff $X^\prime X$ matrices seemed largely to disappear, and references to ridge regression became uncommon, at least in my experience.
In the Statistical Learning course, ridge regression makes a resurgence, with an entirely different purpose. In the era of Big Data, there is a distinct danger that we may overfit data by throwing in lots and lots of predictors (large dimension $m$). One way to mitigate that temptation is to penalize the number or size of the coefficients in the model.
The normal equations minimize the sum of squared deviations in a linear regression model. In other words, the ordinary least squares estimate of $\beta$ solves the problem$$\min_b \Vert Y-X b\Vert^2.$$Given that we are squaring errors, and given that quadratic objective functions have nice smoothness, a tempting way to penalize overfitting is to add a penalty for the size of the coefficient vector:$$\min_b \Vert Y-X b\Vert^2 + \lambda \Vert b \Vert^2$$where $\lambda > 0$ is a penalty parameter. If $\lambda = 0$, we are back to ordinary least squares. As $\lambda \rightarrow \infty$, $\Vert b \Vert$ will converge toward (but never quite reach) 0.
I'll skip the algebra and just state that solving this problem for fixed $\lambda$ is mathematically equivalent to ridge regression as seen above. The only difference is in the choice of $\lambda$. When we were fighting rounding error, we wanted the smallest $\lambda$ that would clean up most of the rounding error. Now, when we are fighting overfitting, we probably need larger values of $\lambda$ to get the job done.
So everything old is new again ... which bodes well for my knees.
During one of the lectures I learned something new that was both interesting and had the side effect of making me feel old (a job my knees are entirely up to on their own, thank you very much). Ridge regression has apparently been picked up from the scrap heap of statistical history and repurposed.
Stiffness: it's not just for us old people
To explain ridge regression as I originally learned it, let me start by defining the phrase ill-conditioned matrix. Let $A$ be a square matrix (we'll operate over the reals, although I think this works as well for complex matrices) and let $\left\Vert \cdot\right\Vert $ be a consistent matrix norm. The condition number of $A$, assuming that is nonsingular, is $$\kappa = \left\Vert A\vphantom{^{-1}}\right\Vert \left\Vert A^{-1}\right\Vert \ge 1.$$ For singular matrices, $\kappa$ is either undefined or $\infty$, according to your religious preferences.
Somewhat informally, $\kappa$ can be interpreted as a measure of how much rounding errors may be inflated then inverting $A$ (using floating-point arithmetic ... "pure" mathematicians don't commit rounding errors). An ill-conditioned or "stiff" matrix $A$ is non-singular but has a large condition number (the criterion for "large" being rather fuzzy), meaning that rounding errors during the process of inverting (or factoring, or otherwise torturing) $A$ may result in an answer with lots of incorrect digits. Long ago, a graduate school instructor told my class the story of an engineer who had inverted an ill-conditioned matrix, without realizing it was ill-conditioned, using home-brew FORTRAN code. He'd published the inverse in a journal article to four or five decimal places. All the decimal places in most of the coefficients were wrong. So were lots of digits to the left of the decimal points.
Ridge regression then
Here comes the tie-in to regression. Suppose we are trying to fit the model $$y=x'\beta+\epsilon$$(where $\epsilon$ is random noise with mean zero), and we have collected observations $(X,Y)$ with $X$ an $n\times m$ matrix ($n$ = sample size, $m$ = number of predictors). The normal equations tell us that the least squares estimate for $\beta$ is $$b=(X^\prime X)^{-1}X^\prime Y.$$In theory, $X^\prime X$ is nonsingular unless the predictors are perfectly multicollinear. In practice, if $X^\prime X$ is nonsingular but stiff, our coefficient estimate $b$ can suffer drastic rounding errors. To mitigate those errors, ridge regression tampers with $X^\prime X$, adding a small positive multiple of the identity matrix of dimension $m$. The modified normal equations are$$b_\lambda=(X^\prime X + \lambda I)^{-1}X^\prime Y.$$The result $b_\lambda$ is a deliberately biased estimate of $\beta$, in which we trade the bias for reduced rounding error. We generally choose $\lambda > 0$ as small as possible subject to the need to stabilize the numerics.
When I first encountered ridge regression, we were in the (first) mainframe era. Machine words were small, single-precision arithmetic ruled the computing universe, and rounding errors were relatively easy to generate. Somewhere along the way, computing hardware got smaller and cheaper, double precision arithmetic became the norm and separate floating point units took over the computational load. At the same time, matrix factorization replaced matrix inversion, and factoring algorithms got smarter about managing rounding error. The problems with stiff $X^\prime X$ matrices seemed largely to disappear, and references to ridge regression became uncommon, at least in my experience.
Ridge regression now
In the Statistical Learning course, ridge regression makes a resurgence, with an entirely different purpose. In the era of Big Data, there is a distinct danger that we may overfit data by throwing in lots and lots of predictors (large dimension $m$). One way to mitigate that temptation is to penalize the number or size of the coefficients in the model.
The normal equations minimize the sum of squared deviations in a linear regression model. In other words, the ordinary least squares estimate of $\beta$ solves the problem$$\min_b \Vert Y-X b\Vert^2.$$Given that we are squaring errors, and given that quadratic objective functions have nice smoothness, a tempting way to penalize overfitting is to add a penalty for the size of the coefficient vector:$$\min_b \Vert Y-X b\Vert^2 + \lambda \Vert b \Vert^2$$where $\lambda > 0$ is a penalty parameter. If $\lambda = 0$, we are back to ordinary least squares. As $\lambda \rightarrow \infty$, $\Vert b \Vert$ will converge toward (but never quite reach) 0.
I'll skip the algebra and just state that solving this problem for fixed $\lambda$ is mathematically equivalent to ridge regression as seen above. The only difference is in the choice of $\lambda$. When we were fighting rounding error, we wanted the smallest $\lambda$ that would clean up most of the rounding error. Now, when we are fighting overfitting, we probably need larger values of $\lambda$ to get the job done.
So everything old is new again ... which bodes well for my knees.
Thursday, August 1, 2013
Numerical Analysis: Size Matters
A pure mathematician, an applied mathematician and an accountant all apply for the same job on Wall Street, and during their interviews they are all asked the same question: "What is the value of $(1-1/N)\times N + 1 - N$?" The pure mathematician replies "Zero, of course!" The applied mathematician asks "How large is $N$?" The accountant asks "What would you like it to be?" The accountant gets the job ...
... but the applied mathematician is on the right track, especially in an era where much research and virtually all application of mathematics is done on digital computers. I've been meaning to do some posts on numerical analysis for a while now, motivated in part by a large number of questions appearing on optimization support forums that indicate the topic is foreign to a sizable number of people using mathematical programming software. (The poster child for questions of this type: "Why did solver XXX say that the optimal value of my variable is 0.999999823 when it is a 0-1 variable??")
I was a math major end to end. In my first quarter of graduate school, I took a course on numerical analysis -- I can't recall now why, as I was a "pure" mathematician in those days -- and it was a real eye-opener. Identities that hold in symbolic mathematics (including "chalkboard arithmetic") don't necessarily survive the conversion to a digital medium. I have no intention of trying to compress that first course in numerical analysis into a blog post, or even a series of posts, but perhaps an example or two will at least get the attention of some people who previously had not heard of it. The numerical analysis course turned out to be one of the most valuable math courses I took (at any level), and I invariably recommend it whenever someone asks "What courses should I take to prepare for a career in Operations Research?"
I tested the formula in the above joke in a small Java program, using double-precision arithmetic, and the reported answer was zero no matter how large a value of $N$ I tried. This might or might not be a consequence of compiler optimizations -- the compiler doing the arithmetic in a different order than I specified it. So I'll fall back on an example from that graduate course: computing a partial sum of the (divergent) harmonic series $\sum_{n=1}^\infty 1/n$.
One of the first things anyone learns about math, at an early age, is that addition iscommutative associative ... except, as it turns out, when a computer is involved. Here is a small Java program that computes $\sum_{n=1}^N 1/n$ twice, in one case summing in ascending order of $n$ ($1 + 1/2 + 1/3 + \dots + 1/n$) and in the other case summing in descending order of $n$ ($1/N + 1/(N-1) + \dots + 1$). I start with $N=1000$ and double $N$ each time through the main loop. To save myself repetition, let me stipulate that henceforth "large" and "small" refer to the magnitudes of numbers, without regard for their sign. If addition is commutative on the computer, the difference between summing in ascending and descending order should be zero.
Here is the output it generated:
The first takeaway is that we're not seeing a difference of zero, ever. (Actually, for $N$ up to about 40, the reported difference is zero.) The differences are small, but not zero. There is no obvious pattern to their sign, nor is their size changing monotonically, but they do seem generally to be growing with $N$, which among other things means that if we have a really large value of $N$, we could get a fairly seriously inaccurate answer (at least in absolute, rather than relative, terms).
Why does this happen? The short answer is that computers store non-integral numbers only to a specific finite precision. That creates "truncation errors" in the representations of the numbers (0.1 as stored in a computer is not exactly 1/10) and "rounding errors" when arithmetic is performed. In particular, adding or subtracting a large number and a small number can produce rounding errors. In the extreme, the small number is smaller than the precision at which the large number is stored, and so the sum or difference of the two is identical to the large number as far as the computer is concerned. That sort of error is more likely to occur when $n$ is ascending (we are adding progressively smaller terms $1/n$ to progressively larger partial sums) than when $n$ is descending (at the outset, we add small numbers to small partial sums, and eventually larger numbers to larger sums); so I am inclined to trust the second set of partial sums more than I do the first. Neither, though, is exactly correct.
You might not be too excited about the differences in the table -- even at $N=1,024,000$ the difference is in the thirteenth decimal place -- but consider the following coding error probably made by every rookie programmer to tackle a mathematical computation: testing whether the difference of two floating point (non-integer) numbers equals zero when they should be testing whether it is close enough to zero. Many a program has entered an indefinite loop thanks to a test for equality and a smidgen of rounding error.
Update: John D. Cook has a related post in a different context (probability and statistics).
... but the applied mathematician is on the right track, especially in an era where much research and virtually all application of mathematics is done on digital computers. I've been meaning to do some posts on numerical analysis for a while now, motivated in part by a large number of questions appearing on optimization support forums that indicate the topic is foreign to a sizable number of people using mathematical programming software. (The poster child for questions of this type: "Why did solver XXX say that the optimal value of my variable is 0.999999823 when it is a 0-1 variable??")
I was a math major end to end. In my first quarter of graduate school, I took a course on numerical analysis -- I can't recall now why, as I was a "pure" mathematician in those days -- and it was a real eye-opener. Identities that hold in symbolic mathematics (including "chalkboard arithmetic") don't necessarily survive the conversion to a digital medium. I have no intention of trying to compress that first course in numerical analysis into a blog post, or even a series of posts, but perhaps an example or two will at least get the attention of some people who previously had not heard of it. The numerical analysis course turned out to be one of the most valuable math courses I took (at any level), and I invariably recommend it whenever someone asks "What courses should I take to prepare for a career in Operations Research?"
I tested the formula in the above joke in a small Java program, using double-precision arithmetic, and the reported answer was zero no matter how large a value of $N$ I tried. This might or might not be a consequence of compiler optimizations -- the compiler doing the arithmetic in a different order than I specified it. So I'll fall back on an example from that graduate course: computing a partial sum of the (divergent) harmonic series $\sum_{n=1}^\infty 1/n$.
One of the first things anyone learns about math, at an early age, is that addition is
for (int limit = 1000; limit <= 1024000; limit *=2) {
double up = 0;
double down = 0;
for (int n = 1; n <= limit; n++) {
up += 1.0/n;
}
for (int n = limit; n >= 1; n--) {
down += 1.0/n;
}
double diff = up - down;
System.out.println("Limit = " + limit + ", up = " + up + ", down = "
+ down + ", diff = " + diff);
}
Here is the output it generated:
| N | 1 + ... + 1/N | 1/N + ... + 1 | Difference |
|---|---|---|---|
| 1000 | 7.485470860550343 | 7.485470860550341 | 2.6645352591003757E-15 |
| 2000 | 8.178368103610284 | 8.178368103610278 | 5.329070518200751E-15 |
| 4000 | 8.871390299795198 | 8.871390299795223 | -2.4868995751603507E-14 |
| 8000 | 9.564474984261391 | 9.564474984261425 | -3.375077994860476E-14 |
| 16000 | 10.257590915797937 | 10.257590915797914 | 2.3092638912203256E-14 |
| 32000 | 10.95072247160203 | 10.950722471602038 | -8.881784197001252E-15 |
| 64000 | 11.64386183972301 | 11.643861839723002 | 8.881784197001252E-15 |
| 128000 | 12.337005114048072 | 12.337005114048194 | -1.2256862191861728E-13 |
| 256000 | 13.030150341487087 | 13.030150341486952 | 1.3500311979441904E-13 |
| 512000 | 13.72329654548546 | 13.723296545485356 | 1.0480505352461478E-13 |
| 1024000 | 14.41644323776362 | 14.416443237764337 | -7.176481631177012E-13 |
The first takeaway is that we're not seeing a difference of zero, ever. (Actually, for $N$ up to about 40, the reported difference is zero.) The differences are small, but not zero. There is no obvious pattern to their sign, nor is their size changing monotonically, but they do seem generally to be growing with $N$, which among other things means that if we have a really large value of $N$, we could get a fairly seriously inaccurate answer (at least in absolute, rather than relative, terms).
Why does this happen? The short answer is that computers store non-integral numbers only to a specific finite precision. That creates "truncation errors" in the representations of the numbers (0.1 as stored in a computer is not exactly 1/10) and "rounding errors" when arithmetic is performed. In particular, adding or subtracting a large number and a small number can produce rounding errors. In the extreme, the small number is smaller than the precision at which the large number is stored, and so the sum or difference of the two is identical to the large number as far as the computer is concerned. That sort of error is more likely to occur when $n$ is ascending (we are adding progressively smaller terms $1/n$ to progressively larger partial sums) than when $n$ is descending (at the outset, we add small numbers to small partial sums, and eventually larger numbers to larger sums); so I am inclined to trust the second set of partial sums more than I do the first. Neither, though, is exactly correct.
You might not be too excited about the differences in the table -- even at $N=1,024,000$ the difference is in the thirteenth decimal place -- but consider the following coding error probably made by every rookie programmer to tackle a mathematical computation: testing whether the difference of two floating point (non-integer) numbers equals zero when they should be testing whether it is close enough to zero. Many a program has entered an indefinite loop thanks to a test for equality and a smidgen of rounding error.
Update: John D. Cook has a related post in a different context (probability and statistics).
Subscribe to:
Posts (Atom)