Kernel Ridge Regression

The ridge page fits a line to the thrown ball and to the crowd by choosing one weight per feature, and it chooses the weights by solving a system as wide as there are features. On the fifteen throws there is one feature, the time, and the line that comes out scores an R² of 0.003, since the ball goes up and comes down and a line can do only one of those.

On the fifteen people the line scores 0.890 with a slope of 0.826 kilograms per centimetre. Nothing on this page changes what ridge minimises.

minimise  ‖y − Xw‖² + λ‖w‖²

w = (XᵀX + λI)⁻¹ Xᵀ y

The matrix being inverted is p by p, where p counts the features. That is the fact to hold onto, because the whole page is about finding a second way to write w whose matrix is n by n, where n counts the rows, and then noticing what the second way lets us do that the first cannot.

Ridge’s answer can be pushed through to a second form. Instead of one weight per feature it becomes one weight per training row, and the fitted w is the sum of the rows, each multiplied by its own weight. The two forms describe the same line, and the derivation is three lines of algebra with one honest condition on it, which section 17 returns to.

w = Xᵀ (XXᵀ + λI)⁻¹ y = Xᵀ a

a = (K + λI)⁻¹ y,   K = XXᵀ

The vector a is called the dual weights, and it has one entry per row. The matrix K is n by n and holds the inner product of every training row with every other, and it is called the Gram matrix. Nothing else about the rows survives into the solve. A prediction is then an inner product too, the new row against each training row, weighted by a.

prediction(x) = xᵀw = Σᵢ aᵢ (x · xᵢ) + ȳ

Take the three people at 169, 170 and 171 centimetres, weighing 64, 70 and 76 kilograms. Centred on their mean of 170 the heights are −1, 0 and 1, and the Gram matrix is the three by three table of their products, which is the only thing the solve ever reads about them. The middle row and column are all zero, because the middle person is at the mean and has no centred height to multiply.

Switch to the radial basis kernel and the same three people produce a different table, with ones down the diagonal, 0.368 for the two pairs a centimetre apart and 0.018 for the pair two centimetres apart. The rows have not changed. What changed is the function that turns two rows into one number, and that function is what the rest of the page varies.

The fifteen throws make a fifteen by fifteen table, and drawn as a heat map the three kernels look nothing alike. The linear kernel’s table is the outer product of one column with itself, bright in the corners and zero along the middle where the centred time passes through nought; the radial one is a band along the diagonal, each throw similar to its neighbours in time and to nothing far away.

The eigenvalues under each table say how many directions the table can tell apart. Under the linear kernel fourteen of the fifteen are zero to within rounding, the largest at about 4e−15, and one is 22.834; a rank-one table can fit exactly one direction, which is a line. The squared polynomial kernel has three nonzero eigenvalues and the radial kernel has fifteen, its smallest 8.1e−9.

A derivation is a claim, and the fifteen throws are where it is checked. Fit ridge regression feature by feature and the dual form row by row under the plain inner product, at the same penalty, and compare everything that can be compared.

At a penalty of one, ridge reports a slope of 0.261135 and the dual weights, added up through the centred times, give 0.261135. Sampled at 81 times across the throw the two lines part by at most 3.7e−14, both score 0.002588, and both answer 12.76 metres at two seconds, which is the mean height because two seconds is the mean time. Drag the penalty through six orders of magnitude and the two columns stay equal to the last printed digit.

The three centred heights are minus one, zero and one; their centred weights are minus six, zero and six. With a ridge penalty of one, the primal calculation is short.

slope = Sxy / (Sxx + λ) = 12 / (2 + 1) = 4

The dual system must recover the same slope through three coefficients, one per training person.

(K + λI) a = y_centred

K + I = [[2, 0, −1], [0, 1, 0], [−1, 0, 2]]

a₂ = 0
2a₁ − a₃ = −6
−a₁ + 2a₃ = 6

solution: a₁ = −2, a₂ = 0, a₃ = 2

(−2)(−1) + (0)(0) + (2)(1) = 4

prediction = mean weight + slope × centred height
           = 70 + 4 × 0.5 = 72 kg

The person at the mean carries a weight of zero and contributes nothing to any prediction, and the two outer people carry equal and opposite weights. A row at the centre of the data tells a straight fit nothing about its slope, and the dual form says so directly where the per-feature form never mentions individual rows at all.

There is a plainer reading of what each weight is, and it holds under every kernel. Rearrange the system that defines the weights and the left side becomes the fitted values, so the weights are the misses of the fit, each divided by the penalty.

K a + λ a = y_centred   ⇒   λ a = y_centred − K a = residuals

aᵢ = residualᵢ / λ

On the three people at a penalty of one the fitted weights are 66, 70 and 74 kilograms, the residuals are −2, 0 and 2, and dividing by one gives the dual weights back. On the fifteen throws the ledger lists every row.

Under the linear kernel at a penalty of one the last two columns agree to 3.2e−14, and under the radial kernel to 2.7e−15. Drag the penalty down and something worth noticing happens under the linear kernel. The line barely moves, since its score is 0.002593 at a thousandth and 0.002588 at one, while the largest dual weight goes from 12.08 to 12,105, growing as one over the penalty because the residuals stay where they are and the divisor shrinks.

Ridge regression on the crowd keeps three numbers, a weight for height, an intercept, and nothing else, and “0.826 kilograms per centimetre” is a sentence about height. The dual model keeps every training row and a weight for each, so on fifteen people it keeps fifteen heights and fifteen weights, and a prediction is a sum over all of them.

There is no per-feature number to report. It is possible to add the weights up through the centred heights and recover the slope, as section 5 did, only because under the linear kernel there is a slope to recover; under the radial kernel there is not, and the playground says so rather than inventing one.

the model = the n training rows + the n dual weights + ȳ

Everything in Part 1 reads the rows through one function, the inner product of two of them. The kernel trick page argues that any function which is an inner product in some expanded space, whether or not that space is ever built, can stand in for it, and this page takes that as given. Replace the inner product wherever it appears, in the Gram matrix and in the prediction, and nothing else about the fit changes.

K[i, j] = k(xᵢ, xⱼ)

a = (K + λI)⁻¹ y_centred

prediction(x) = Σᵢ aᵢ k(x, xᵢ) + ȳ

The solve is the same solve, the penalty is the same penalty on the same diagonal, and the dual weights are still the residuals over the penalty. What the swap buys is a curve, drawn without a single manufactured column, in the original picture of time against height.

The polynomial kernel raises the inner product plus one to a power, and the plus one is what lets every lower power ride along with the highest. On the fifteen throws the degree does what the polynomial page’s manufactured columns did, except that no column is ever built.

…

Degree 2 follows the arc at 0.981, degree 3 reaches 0.985, and degrees 4, 5 and 6 come out at 0.981, 0.982 and 0.980, so past the square the extra powers buy nothing on this throw at this penalty. What they cost is visible in the condition number of the system, which goes 23.8, 72.4, 262, 1061, 4597 and 20,781 as the degree rises, since the largest entries of the Gram matrix are centred times raised to the degree and the smallest stay near zero.

The textbook reading is that the polynomial kernel at degree one is not the linear kernel, because the plus one adds a constant to every entry of the Gram matrix and a constant entry is a feature that is always one. I measured it and found the two fits identical here. On the throw at a penalty of one the two curves part by 2.0e−14, and on the three people by 3.9e−12, which is the arithmetic of a solve and not a difference of model.

K_degree1 = K_linear + 11ᵀ

The radial basis kernel reads two rows as similar when they are close, one at zero distance and falling toward zero as the distance grows, with gamma setting how fast. The row of the Gram matrix for the throw at two seconds reads 1.0 at itself, then 0.919, 0.723, 0.477, 0.273, 0.129, 0.054 and 0.018 for the throws stepping outward at gamma 1, so each throw is compared mainly with its neighbours in time.

…

At gamma 0.01 every throw looks like every other and the fit is a curve that is nearly straight at 0.012. At 0.1 it reaches 0.450, at 0.3 it follows the arc at 0.797, at 1 it scores 0.873, at 3 it peaks at 0.880, and by 10 each throw is similar only to itself, the smallest Gram eigenvalue has risen from 8.1e−9 at gamma 1 to 0.202, and the curve dips toward the mean between the throws for a score of 0.844.

At a hundred the table is the identity to four decimals and the fit is the mean plus a bump at each throw.

Gamma multiplies a squared distance, so it carries the inverse square of whatever the feature is measured in. The gamma of 1 that suited times in seconds is a reach of one second; applied to heights in centimetres it is a reach of one centimetre, and the people are three apart.

…

At gamma 1 on the crowd the Gram matrix is the identity to four decimals, the smallest eigenvalue 0.99976, so every dual weight is half that person’s centred weight and the curve visits each person and returns to the mean of 72 between them. Asked about the person at 173 centimetres it answers 74.0, halfway from the mean to their 76 kilograms, and asked about 174.5, between two people, it answers 72.16.

The score is 0.750, below the line’s 0.890. Widen the reach to a hundredth and the score is 0.835; at a thousandth 0.824; at a ten-thousandth the kernel can no longer tell the tallest from the shortest and the score falls to 0.481. Nothing here beats the line, because the crowd was drawn along one.

The three people are small enough that every kernel’s answer can be read whole. Each fit is asked the same question, what a person of 170.5 centimetres weighs, and each answers from its own three dual weights.

The middle weight is zero under every kernel, for the reason section 6 gave. The middle person weighs the mean, every one of these fits passes through the mean, so their residual is zero and so is their weight. The outer weights differ because the kernels differ about how much the outer two people resemble each other, 1 against −1 under the linear kernel and 0.018 under the radial one.

On the ridge page the penalty matrix has a zero in its first slot so the intercept is not shrunk, and the page records the bug that came from zeroing that slot when there was no intercept column to exempt. The dual system has no columns at all. Its matrix is indexed by training rows on both sides, and no row is an intercept, so the penalty is added to every entry of the diagonal with no exemption.

K + λI,   every diagonal entry, no exemption

Exempting row zero would leave one arbitrary person unpenalised, a fit that still runs, still draws a curve, and is quietly wrong. What does the intercept’s job instead is the centring. Both the rows and the target have their means removed before the solve and the target mean is added back to every prediction, which is why every fit on this page passes through the mean point and why the throw answers exactly 12.76 at two seconds under any kernel.

The ridge page describes the penalty as a trade of fit against restraint. In the dual there is a second thing it does that has no counterpart on that page, since the Gram matrix of a radial kernel on distinct rows is positive definite and only barely, its smallest eigenvalue on the throw at gamma 1 being 8.1e−9, and adding the penalty to the diagonal is what makes the solve possible at all.

…

At a penalty of a billionth the radial curve scores 0.999998 on the throws and swings 122 metres from the mean between them, the largest dual weight is 1.6e7 and the condition number of the system is 6.2e8, which is nine digits of accuracy lost from sixteen. At a millionth the score is still 0.99955 with a condition number of 5.6e6 and a swing of 17.8; at a tenth it is 0.987 with a condition number of 57.5 and the curve stays within 11 metres of the mean; at one, 0.873 and 6.65.

condition number = (largest eigenvalue + λ) / (smallest eigenvalue + λ)

at λ = 1:    condition number ≈ 6.650 / 1 ≈ 6.65
at λ = 10⁻⁹: condition number ≈ 5.650 / (8.1 × 10⁻⁹ + 10⁻⁹)
                             ≈ 6.2 × 10⁸

At the other end the penalty dominates and the curve goes flat. With the dual weights being residuals over the penalty, a large penalty makes every weight small whatever the residuals are, and a prediction is the mean plus a sum of small numbers.

By a thousand the curve is within three centimetres of the mean everywhere across the throw and the score is 0.005. Under the linear kernel the same thing happens more slowly, its score going from 0.002588 at one to 0.000114 at a thousand, because there was almost nothing to flatten.

The derivation in section 2 holds on centred data. Both the target and the rows have to have their means removed, and when this model was first built here only the target was. The result is the most instructive kind of bug, a fit that still ran, still drew a plausible line and agreed with ridge to within two percent, 4.2841 against 4.3580 on the fixture it was first measured on.

On the crowd the mistake is not subtle at all. The heights are a hundred and seventy centimetres from the origin, so the uncentred Gram matrix is dominated by that offset, and the half-centred fit reports a slope of 0.0046 kilograms per centimetre where both right fits report 0.8259, and answers 72.80 at the mean height where the right lines answer exactly 72.

Drop the penalty to a thousandth and the right fits rise to 0.8262 while the wrong one stays at 0.0046. On the throw, where the times are near the origin, the wrong slope is 0.0742 against 0.2611, which is the closer-looking version of the same failure.

The trade the two forms set up is a matrix as wide as the features against a matrix as wide as the rows, and it is a claim about time, so it was timed. Both fits run on drawn data, first with the rows doubling at two features and then with the features doubling at two hundred rows, each timing the median of five runs through the whole fit, boundary checks included.

… timing the fits, which takes a moment the first time

On the machine I measured on, the kernel fit took 0.09 milliseconds at fifty rows and 121 milliseconds at sixteen hundred, roughly a thousandfold for a thirty-twofold increase in rows, which is the n-cubed solve of an n by n system, while ridge stayed between 0.02 and 0.06 milliseconds throughout. Growing the features instead, from two to sixty-four at two hundred rows, took ridge from 0.04 to 0.30 milliseconds and the kernel fit from 1.2 to 3.1, the kernel’s share coming from building the Gram matrix, which touches every feature once per pair of rows.

ridge:  solve p × p,   kernel ridge:  solve n × n

Prediction has the same shape. A ridge prediction multiplies p weights; a kernel prediction computes the kernel between the query and every training row and multiplies n weights, so predicting a hundred rows took 0.095 milliseconds against fifty training rows and 2.85 against sixteen hundred, where ridge stayed near 0.015. What the fitted model has to keep grows the same way, the rows times the features plus one weight per row, 150 numbers at fifty rows and 4800 at sixteen hundred, against ridge’s three.

A radial kernel is a bump around each training row, and a query far from every row is similar to none of them, so every term in its prediction is near zero and what is left is the mean. The throw was measured for four seconds. Ask the radial fit what the ball is doing at six.

At gamma 1 and a penalty of a tenth the curve answers 7.47 metres at five seconds, already off the last throw, 12.46 at six, 12.758 at seven and 12.759998 at eight, against a mean of 12.76. Two seconds past the data the fit has forgotten the ball entirely and is answering the average height of the throw.

The squared polynomial kernel does the opposite and keeps its shape, answering −19.5 at five seconds, −49.8 at six and −136.5 at eight, which is the polynomial page’s runaway drawn without columns. The line answers 14.33 at eight seconds and would keep the same slope forever.

Two things ridge reports that this model cannot. A coefficient per feature, for the reason section 7 gave, and an intercept, because there was never an intercept column. What there is instead is the target mean, added back to every prediction, and on the crowd it is 72 kilograms. Under the linear kernel a slope and an intercept can be reconstructed from the weights and the mean, and the centring widget does exactly that to compare against ridge; under any other kernel the reconstruction has nothing to reconstruct.

Not every function of two rows is a kernel. A kernel’s Gram matrix has no negative eigenvalue on any set of rows, which is Mercer’s condition, and the sigmoid kernel fails it for many settings. The solve uses a Cholesky factorisation, which succeeds exactly when the matrix it is given is positive definite, so the factorisation is also the test, and the fit is refused at the first pivot that is not positive.

On the throws at gamma 1 and a constant of −0.5 the smallest eigenvalue of the Gram matrix is −4.313. At every penalty up to 4.3 the sum with the smallest eigenvalue is negative and the fit is refused, in a message naming Mercer’s condition. At 4.4 the sum is positive, the factorisation goes through, and the fit is accepted with a score of −228, which is a curve far worse than answering the mean; at 5 it scores −4.49, at 10 −0.069, and at 100 it is 0.006, the flat line the penalty forces.

Even at a constant of zero the smallest eigenvalue is −0.878, hidden at a penalty of one, and the accepted fit scores 0.0023.

Under the linear kernel on the fifteen throws the Gram matrix has one nonzero eigenvalue and fourteen at zero, because fifteen rows of one centred feature span one direction. Without the penalty the system to be solved is that rank-one table, which has no unique solution, and under the squared polynomial kernel twelve of the fifteen eigenvalues are zero for the same reason. The radial kernel’s table has no exact zero, its smallest at 8.1e−9, and section 15 showed what solving near it costs.

λ = 0:  K a = y_centred has no unique a when K is singular

A penalty of zero is refused when the model is constructed, before any data is seen, and this page refuses it a layer earlier still. The sweep in section 15 stops at a billionth for the same reason, since that is where the solve still completes and the condition number is already 6.2e8.

A complete implementation states which kernel and which of its parameters, whether it centres the rows and the target and where the target mean is put back, that the penalty is added to every diagonal entry, which solver factorises the system and what it does when the factorisation fails, whether the dual weights are exposed and that they are not coefficients, that the training rows are kept and are part of the saved model, how the query features are matched to the fitted ones, and the smallest penalty it will accept.

The model here centres both, adds the penalty everywhere, solves by Cholesky, keeps the rows, matches features by name, exposes the dual weights and the training rows as copies, and refuses a penalty that is not positive at construction. Every row of the table below was run against it.

The constant column deserves the extra sentence, because the two models on this page treat it oppositely. The kernel model alone accepts it and answers the mean everywhere, since with no spread in the column there is nothing for a kernel to compare; the ridge control fitted beside it for the comparison refuses, so the page’s endpoint refuses, and the refusal belongs to ridge rather than to the dual.