linear regression with parameters constrained to lay on the unit ball

lr.R

# fake data.
X = matrix(data = rnorm(20), nr = 5, nc = 4);
y = X %*% (runif(4) * 2 - 1) + rnorm(5) * 0.1

# the usual linear regression.
b = solve(t(X) %*% X) %*% t(X) %*% y
sum((y - X %*% b)^2)
sum(b*b)

# the unit ball constrained linear regression.
v = t(X) %*% y
E = eigen(t(X) %*% X)
uv = t(E$vectors) %*% v
lam_min = -E$values[4];
lam_max = sum(uv^2);
while(lam_max - lam_min > 0.001) {
	lam_mid = (lam_max + lam_min) / 2
	f = sum(uv^2 * (E$values+lam_mid)^(-2)) - 1
	if(f < 0) {
		lam_max = lam_mid
	} else {
		lam_min = lam_mid
	}
	print(c(f, lam_min, lam_max))
}
lam = (lam_max + lam_min) / 2
c = solve(t(X) %*% X + lam*diag(4)) %*% v
sum((y - X %*% c)^2)
sum(c*c)

# the projected thing
d = b / sqrt(sum(b*b))
sum((y - X %*% d)^2)
sum(d*d)