# How to calculate regularization parameter in ridge regression given degrees of freedom and input matrix?

Let A be $n \times p$ matrix of independent variables and B be the corresponding $n \times 1$ matrix of the dependent values. In ridge regression, we define a parameter $\lambda$ so that: $\beta=(A^\mathrm{T}A+\lambda I)^{-1}A^\mathrm{T}B$ . Now let [u s v]=svd(A) and $d_{i}=i^{th}$ diagonal entry of 's'. we define degrees of freedom (df)= $\sum_{i=1}^{n} \frac{(d_{i})^2}{(d_{i})^2+\lambda}$ . Ridge regression shrinks the coefficients of low-variance components and hence the parameter $\lambda$ controls the degrees of freedom.So for $\lambda=0$, which is the case of normal regression, df=n, and hence all the independent variables will be considered. The problem I am facing is to find the value of $\lambda$ given 'df' and the matrix 's'. I have tried to re-arrange the above equation but was not getting a closed form solution. Please provide any helpful pointers.

• Well I need time to answer this (probably others will be quicker to aid you), but most insights may be taken from stat.lsa.umich.edu/~kshedden/Courses/Stat600/Notes/… And what is $k$ in definition of degrees of freedom, since I do miss $\lambda$ somehow. – Dmitrij Celov Mar 15 '11 at 10:55
• @Dmitrij: Thnx for the reply, I have updated the questions, and replaced 'k' with $\lambda$ – Amit Mar 15 '11 at 11:23
• Hi Amit, how can you know what the degrees of freedom are before calculating the regularization parameter? – Baz Oct 11 '13 at 16:22

A Newton-Raphson/Fisher-scoring/Taylor-series algorithm would be suited to this.

You have the equation to solve for $\lambda$ $$h(\lambda)=\sum_{i=1}^{p}\frac{d_{i}^{2}}{d_{i}^{2}+\lambda}-df=0$$ with derivative $$\frac{\partial h}{\partial \lambda}=-\sum_{i=1}^{p}\frac{d_{i}^{2}}{(d_{i}^{2}+\lambda)^{2}}$$ You then get: $$h(\lambda)\approx h(\lambda^{(0)})+(\lambda-\lambda^{(0)})\frac{\partial h}{\partial \lambda}|_{\lambda=\lambda^{(0)}}=0$$

re-arranging for $\lambda$ you get: $$\lambda=\lambda^{(0)}-\left[\frac{\partial h}{\partial \lambda}|_{\lambda=\lambda^{(0)}}\right]^{-1}h(\lambda^{(0)})$$ This sets up the iterative search. For initial starting values, assume $d^{2}_{i}=1$ in the summation, then you get $\lambda^{(0)}=\frac{p-df}{df}$.

$$\lambda^{(j+1)}=\lambda^{(j)}+\left[\sum_{i=1}^{p}\frac{d_{i}^{2}}{(d_{i}^{2}+\lambda^{(j)})^{2}}\right]^{-1}\left[\sum_{i=1}^{p}\frac{d_{i}^{2}}{d_{i}^{2}+\lambda^{(j)}}-df\right]$$

This "goes" in the right direction (increase $\lambda$ when summation is too big, decrease it when too small), and typically only takes a few iterations to solve. Further the function is monotonic (an increase/decrease in $\lambda$ will always decrease/increase the summation), so that it will converge uniquely (no local maxima).

• thnx a lot , but I have a doubt why do we need to assume $d_{i}^2=1$, since we have their correct values already...I have checked this formula by writing a matlab code and have not taken that assumption, but it works fine and gives correct solution – Amit Mar 15 '11 at 13:21
• The assumption is just to get the initial value of $\lambda^{(0)}$ "close enough" to the correct value. If you have a better guess, then start with that. You could even just set $\lambda^{(0)}=0$, as long as your d's are greater than zero. The d's are not assumed 1 in the iterations, just to get the algorithm started. – probabilityislogic Mar 15 '11 at 13:27
• (+1) I would give the same numerical solution anyway. – Dmitrij Celov Mar 15 '11 at 14:44

Here is the small Matlab code based on the formula proved by probabilityislogic:

function [lamda] = calculate_labda(Xnormalised,df)
[n,p] = size(Xnormalised);

%Finding SVD of data
[u s v]=svd(Xnormalised);
Di=diag(s);
Dsq=Di.^2;

%Newton-rapson method to solve for lamda
lamdaPrev=(p-df)/df;
lamdaCur=Inf;%random large value
diff=lamdaCur-lamdaPrev;
threshold=eps(class(XstdArray));
while (diff>threshold)
numerator=(sum(Dsq ./ (Dsq+lamdaPrev))-df);
denominator=sum(Dsq./((Dsq+lamdaPrev).^2));
lamdaCur=lamdaPrev+(numerator/denominator);
diff=lamdaCur-lamdaPrev;
lamdaPrev=lamdaCur;
end
lamda=lamdaCur;
end

• Go team! – probabilityislogic Mar 17 '11 at 13:21
• An attempted editor argues that the while condition should be while ( abs(diff)>threshold ). – gung Oct 17 '18 at 13:24
• I am posting this as an alternative answer to the code posted by @Amit. I am suggesting that comparison to threshold in the while( abs(diff) > threshold ) because the tolerance for the difference should be reachable from both the left and the right. For example let's say diff = $-100$ and threshold = $1e^{-16}$ then the loop condition is false, but clearly the solution has not converged. – heuristicwondering Oct 28 '18 at 23:01