405 lines
12 KiB
Julia
405 lines
12 KiB
Julia
|
|
|
|
|
|
"""
|
|
Fast k-fold cv for updating regression coefficients
|
|
"""
|
|
function TRSegCVUpdate(X, y, lambdas, cv, bOld, regType="L2", derOrder=0)
|
|
|
|
n, p = size(X);
|
|
|
|
# Finding appropriate regularisation matrix
|
|
if regType == "bc"
|
|
regMat = [I(p); zeros(derOrder,p)];
|
|
for i = 1:derOrder regMat = diff(regMat, dims = 1); end
|
|
elseif regType == "legendre"
|
|
regMat = [I(p); zeros(derOrder,p)];
|
|
for i = 1:derOrder regMat = diff(regMat, dims = 1); end
|
|
P, _ = plegendre(derOrder-1, p);
|
|
epsilon = 1e-14;
|
|
regMat[end-derOrder+1:end,:] = sqrt(epsilon) * P;
|
|
elseif regType == "L2"
|
|
regMat = I(p);
|
|
elseif regType == "std"
|
|
regMat = Diagonal(vec(std(X, dims=1)));
|
|
elseif regType == "GL" # GL fractional derivative regulariztion
|
|
C = ones(p)*1.0;
|
|
for k in 2:p
|
|
C[k] = (1-(derOrder+1)/(k-1)) * C[k-1];
|
|
end
|
|
|
|
regMat = zeros(p,p);
|
|
|
|
for i in 1:p
|
|
regMat[i:end, i] = C[1:end-i+1];
|
|
end
|
|
end
|
|
|
|
# Preliminary calculations
|
|
mX = mean(X, dims=1);
|
|
X = X .- mX;
|
|
my = mean(y);
|
|
y = y .- my;
|
|
X = X / regMat;
|
|
U, s, V = svd(X, full=false);
|
|
n_seg = maximum(cv);
|
|
n_lambdas = length(lambdas);
|
|
|
|
# Finding residuals
|
|
denom = broadcast(.+, broadcast(./, lambdas, s'), s')';
|
|
denom2 = broadcast(.+, ones(n), broadcast(./, lambdas', s.^2))
|
|
yhat = broadcast(./, s .* (U'*y), denom)
|
|
yhat += s .* broadcast(.-, 1, broadcast(./, 1, denom2)) .* (V' * bOld)
|
|
resid = broadcast(.-, y, U * yhat)
|
|
|
|
# Finding cross-validated residuals
|
|
rescv = zeros(n, n_lambdas);
|
|
sdenom = sqrt.(broadcast(./, s, denom))';
|
|
|
|
for seg in 1:n_seg
|
|
|
|
Useg = U[vec(cv .== seg),:];
|
|
Id = 1.0 * I(size(Useg,1)) .- 1/n;
|
|
|
|
for k in 1:n_lambdas
|
|
Uk = Useg .* sdenom[k,:]';
|
|
rescv[vec(cv .== seg), k] = (Id - Uk * Uk') \ resid[vec(cv .== seg), k];
|
|
end
|
|
end
|
|
|
|
# Calculating rmsecv and regression coefficients
|
|
press = sum(rescv.^2, dims=1)';
|
|
rmsecv = sqrt.(1/n .* press);
|
|
bcoeffs = V * broadcast(./, (U' * y), denom);
|
|
bcoeffs = regMat \ bcoeffs;
|
|
|
|
# Creating regression coefficients for uncentred data
|
|
if my != 0
|
|
bcoeffs = [my .- mX*bcoeffs; bcoeffs];
|
|
end
|
|
|
|
# Finding rmsecv-minimal lambda value and associated regression coefficients
|
|
lambda_min, lambda_min_ind = findmin(rmsecv)
|
|
lambda_min_ind = lambda_min_ind[1]
|
|
b_lambda_min = bcoeffs[:,lambda_min_ind]
|
|
|
|
return b_lambda_min, rmsecv, lambda_min, lambda_min_ind
|
|
end
|
|
|
|
"""
|
|
Updates regression coefficient by solving the augmented TR problem [Xs; sqrt(lambda)*I] * beta = [ys; sqrt(lambda)*b_old]
|
|
Note that many regularization types are supported but the regularization is on the difference between new and old reg. coeffs.
|
|
and so most regularization types are probably not meaningful.
|
|
"""
|
|
function TRLooCVUpdate(X, y, lambdas, bOld, regType="L2", derOrder=0)
|
|
|
|
n, p = size(X);
|
|
|
|
# Finding appropriate regularisation matrix
|
|
if regType == "bc"
|
|
regMat = [I(p); zeros(derOrder,p)];
|
|
for i = 1:derOrder regMat = diff(regMat, dims = 1); end
|
|
elseif regType == "legendre"
|
|
regMat = [I(p); zeros(derOrder,p)];
|
|
for i = 1:derOrder regMat = diff(regMat, dims = 1); end
|
|
P, _ = plegendre(derOrder-1, p);
|
|
epsilon = 1e-14;
|
|
regMat[end-derOrder+1:end,:] = sqrt(epsilon) * P;
|
|
elseif regType == "L2"
|
|
regMat = I(p);
|
|
elseif regType == "std"
|
|
regMat = Diagonal(vec(std(X, dims=1)));
|
|
elseif regType == "GL" # GL fractional derivative regulariztion
|
|
C = ones(p)*1.0;
|
|
for k in 2:p
|
|
C[k] = (1-(derOrder+1)/(k-1)) * C[k-1];
|
|
end
|
|
regMat = zeros(p,p);
|
|
for i in 1:p
|
|
regMat[i:end, i] = C[1:end-i+1];
|
|
end
|
|
end
|
|
|
|
|
|
# Preliminary calculations
|
|
mX = mean(X, dims=1);
|
|
X = X .- mX;
|
|
my = mean(y);
|
|
y = y .- my;
|
|
X = X / regMat;
|
|
U, s, V = svd(X, full=false);
|
|
n_seg = maximum(cv);
|
|
n_lambdas = length(lambdas);
|
|
|
|
# Main calculations
|
|
denom = broadcast(.+, broadcast(./, lambdas, s'), s')';
|
|
denom2 = broadcast(.+, ones(n), broadcast(./, lambdas', s.^2))
|
|
yhat = broadcast(./, s .* (U'*y), denom);
|
|
yhat += s .* broadcast(.-, 1, broadcast(./, 1, denom2)) .* (V' * bOld)
|
|
resid = broadcast(.-, y, U * yhat)
|
|
H = broadcast(.+, U.^2 * broadcast(./, s, denom), 1/n);
|
|
rescv = broadcast(./, resid, broadcast(.-, 1, H));
|
|
press = vec(sum(rescv.^2, dims=1));
|
|
rmsecv = sqrt.(1/n .* press);
|
|
gcv = vec(broadcast(./, sum(resid.^2, dims=1), mean(broadcast(.-, 1, H), dims=1).^2));
|
|
|
|
# Finding lambda that minimises rmsecv and GCV
|
|
idminrmsecv = findmin(press)[2][1];
|
|
idmingcv = findmin(gcv)[2][1];
|
|
lambdarmsecv = lambdas[idminrmsecv];
|
|
lambdagcv = lambdas[idmingcv];
|
|
|
|
# Calculating regression coefficients
|
|
bcoeffs = V * broadcast(./, (U' * y), denom);
|
|
bcoeffs = regMat \ bcoeffs;
|
|
|
|
if my != 0
|
|
bcoeffs = [my .- mX*bcoeffs; bcoeffs];
|
|
end
|
|
|
|
brmsecv = bcoeffs[:, idminrmsecv];
|
|
bgcv = bcoeffs[:, idmingcv];
|
|
|
|
|
|
return brmsecv, bgcv, rmsecv, gcv, idminrmsecv, lambdarmsecv, idmingcv, lambdagcv
|
|
end
|
|
|
|
|
|
"""
|
|
### TO DO: ADD FRACTIONAL DERIVATIVE REGULARIZATION <-- Check that it is correctly added:) ###
|
|
|
|
regularizationMatrix(X; regType="L2", regParam1=0, regParam2=1e-14)
|
|
regularizationMatrix(p::Int64; regType="L2", regParam1=0, regParam2=1e-14)
|
|
|
|
Calculates and returns square regularization matrix.
|
|
|
|
Inputs:
|
|
- X/p : Size of regularization matrix or data matrix (size of reg. mat. will then be size(X,2)
|
|
- regType : "L2" (returns identity matrix)
|
|
"bc" (boundary condition, forces zero on right endpoint for derivative regularization) or
|
|
"legendre" (no boundary condition, but fills out reg. mat. with lower order polynomial trends to get square matrix)
|
|
"std" (standardization, FILL OUT WHEN DONE)
|
|
- regParam1 : Int64, Indicates degree of derivative regularization (0 gives L\\_2)
|
|
- regParam2 : For regType=="plegendre" added polynomials are multiplied by sqrt(regParam2)
|
|
|
|
Output: Square regularization matrix
|
|
"""
|
|
function regularizationMatrix(X; regType="L2", regParam1=0, regParam2=1e-14)
|
|
|
|
if regType == "std"
|
|
regParam2 = Diagonal(vec(std(X, dims=1)));
|
|
end
|
|
|
|
regMat = regularizationMatrix(size(X,2); regType, regParam1, regParam2)
|
|
|
|
return regMat
|
|
end
|
|
|
|
function regularizationMatrix(p::Int64; regType="L2", regParam1=0, regParam2=1e-14)
|
|
|
|
if regType == "bc" # Discrete derivative with boundary conditions
|
|
regMat = [I(p); zeros(regParam1,p)]; for i = 1:regParam1 regMat = diff(regMat, dims = 1); end
|
|
elseif regType == "legendre" # Fill in polynomials in bottom row(s) to get square matrix
|
|
regMat = [I(p); zeros(regParam1,p)]; for i = 1:regParam1 regMat = diff(regMat, dims = 1); end
|
|
P, _ = plegendre(regParam1-1, p);
|
|
regMat[end-regParam1+1:end,:] = sqrt(regParam2) * P;
|
|
elseif regType == "L2"
|
|
regMat = I(p);
|
|
elseif regType == "std" # Standardization
|
|
regMat = regParam2;
|
|
elseif regType == "GL" # Grünwald-Letnikov fractional derivative regulariztion
|
|
# regParam1 is alpha (order of fractional derivative)
|
|
C = ones(p)*1.0;
|
|
for k in 2:p
|
|
C[k] = (1-(regParam1+1)/(k-1)) * C[k-1];
|
|
end
|
|
|
|
regMat = zeros(p,p);
|
|
|
|
for i in 1:p
|
|
regMat[i:end, i] = C[1:end-i+1];
|
|
end
|
|
end
|
|
|
|
return regMat
|
|
end
|
|
|
|
"""
|
|
function TRLooCV
|
|
|
|
bpress, bgcv, rmsecv, GCV, idminPRESS, lambdaPRESS, idminGCV, lambdaGCV = TRLooCV(X, y, lambdas, regType="L2", regParam1=1, regParam2=1)
|
|
|
|
regType: 'bc', 'legendre', 'L2', 'std', 'GL'
|
|
"""
|
|
function TRLooCV(X, y, lambdas, regType="L2", regParam1=1, regParam2=1)
|
|
|
|
n, p = size(X);
|
|
mX = mean(X, dims=1);
|
|
X = X .- mX;
|
|
my = mean(y);
|
|
y = y .- my;
|
|
|
|
if regType == "bc"
|
|
regMat = [I(p); zeros(regParam1,p)]; for i = 1:regParam1 regMat = diff(regMat, dims = 1); end
|
|
elseif regType == "legendre"
|
|
regMat = [I(p); zeros(regParam1,p)]; for i = 1:regParam1 regMat = diff(regMat, dims = 1); end
|
|
P, _ = plegendre(regParam-1, p);
|
|
regMat[end-regParam1+1:end,:] = sqrt(regParam2) * P;
|
|
elseif regType == "L2"
|
|
regMat = I(p);
|
|
elseif regType == "std"
|
|
regMat = Diagonal(vec(std(X, dims=1)));
|
|
elseif regType == "GL" # Grünwald-Letnikov fractional derivative regulariztion
|
|
# regParam1 is alpha (order of fractional derivative)
|
|
C = ones(p)*1.0;
|
|
for k in 2:p
|
|
C[k] = (1-(regParam1+1)/(k-1)) * C[k-1];
|
|
end
|
|
|
|
regMat = zeros(p,p);
|
|
|
|
for i in 1:p
|
|
regMat[i:end, i] = C[1:end-i+1];
|
|
end
|
|
end
|
|
|
|
|
|
X = X / regMat;
|
|
U, s, V = svd(X, full=false)
|
|
|
|
denom = broadcast(.+, broadcast(./, lambdas, s'), s')';
|
|
H = broadcast(.+, U.^2 * broadcast(./, s, denom), 1/n);
|
|
resid = broadcast(.-, y, U * broadcast(./, s .* (U'*y), denom));
|
|
rescv = broadcast(./, resid, broadcast(.-, 1, H));
|
|
press = vec(sum(rescv.^2, dims=1));
|
|
rmsecv = sqrt.(1/n .* press);
|
|
GCV = vec(broadcast(./, sum(resid.^2, dims=1), mean(broadcast(.-, 1, H), dims=1).^2));
|
|
|
|
idminPRESS = findmin(press)[2][1]; # First index selects coordinates, second selects '1st coordinate'
|
|
idminGCV = findmin(GCV)[2][1]; # First index selects coordinates, second selects '1st coordinate'
|
|
lambdaPRESS = lambdas[idminPRESS];
|
|
lambdaGCV = lambdas[idminGCV];
|
|
|
|
bcoeffs = V * broadcast(./, (U' * y), denom);
|
|
bcoeffs = regMat \ bcoeffs;
|
|
|
|
if my != 0
|
|
bcoeffs = [my .- mX*bcoeffs; bcoeffs];
|
|
end
|
|
|
|
bpress = bcoeffs[:, idminPRESS]
|
|
bgcv = bcoeffs[:, idminGCV]
|
|
|
|
return bpress, bgcv, rmsecv, GCV, idminPRESS, lambdaPRESS, idminGCV, lambdaGCV
|
|
end
|
|
|
|
"""
|
|
function TRSegCV(X, y, lambdas, cv, regType="L2", regParam1=0, regParam2=1e-14)
|
|
|
|
Segmented cross-validation based on the Sherman-Morrison-Woodbury updating formula.
|
|
Inputs:
|
|
- X : Data matrix
|
|
- y : Response vector
|
|
- lambdas : Vector of regularization parameter values
|
|
- cv : Vector of length n indicating segment membership for each sample
|
|
- regType, regParam1, regParam2 : Inputs to regularizationMatrix function
|
|
|
|
Outputs: b_lambda_min, rmsecv, lambda_min, lambda_min_ind
|
|
"""
|
|
function TRSegCV(X, y, lambdas, cv, regType="L2", regParam1=0, regParam2=1e-14)
|
|
|
|
n, p = size(X);
|
|
mX = mean(X, dims=1);
|
|
X = X .- mX;
|
|
my = mean(y);
|
|
y = y .- my;
|
|
|
|
if regType == "bc"
|
|
regMat = [I(p); zeros(regParam1,p)]; for i = 1:regParam1 regMat = diff(regMat, dims = 1); end
|
|
elseif regType == "legendre"
|
|
regMat = [I(p); zeros(regParam1,p)]; for i = 1:regParam1 regMat = diff(regMat, dims = 1); end
|
|
P, _ = plegendre(regParam-1, p);
|
|
regMat[end-regParam1+1:end,:] = sqrt(regParam2) * P;
|
|
elseif regType == "L2"
|
|
regMat = I(p);
|
|
elseif regType == "std"
|
|
regMat = Diagonal(vec(std(X, dims=1)));
|
|
elseif regType == "GL" # Grünwald-Letnikov fractional derivative regulariztion
|
|
# regParam1 is alpha (order of fractional derivative)
|
|
C = ones(p)*1.0;
|
|
for k in 2:p
|
|
C[k] = (1-(regParam1+1)/(k-1)) * C[k-1];
|
|
end
|
|
|
|
regMat = zeros(p,p);
|
|
|
|
for i in 1:p
|
|
regMat[i:end, i] = C[1:end-i+1];
|
|
end
|
|
end
|
|
|
|
|
|
X = X / regMat;
|
|
U, s, V = svd(X, full=false);
|
|
|
|
n_seg = maximum(cv);
|
|
n_lambdas = length(lambdas);
|
|
|
|
denom = broadcast(.+, broadcast(./, lambdas, s'), s')';
|
|
resid = broadcast(.-, y, U * broadcast(./, s .* (U'*y), denom));
|
|
rescv = zeros(n, n_lambdas);
|
|
sdenom = sqrt.(broadcast(./, s, denom))';
|
|
|
|
for seg in 1:n_seg
|
|
|
|
Useg = U[vec(cv .== seg),:];
|
|
Id = 1.0 * I(size(Useg,1)) .- 1/n;
|
|
|
|
for k in 1:n_lambdas
|
|
Uk = Useg .* sdenom[k,:]';
|
|
rescv[vec(cv .== seg), k] = (Id - Uk * Uk') \ resid[vec(cv .== seg), k];
|
|
end
|
|
end
|
|
|
|
press = sum(rescv.^2, dims=1)';
|
|
rmsecv = sqrt.(1/n .* press);
|
|
bcoeffs = V * broadcast(./, (U' * y), denom);
|
|
bcoeffs = regMat \ bcoeffs;
|
|
|
|
if my != 0
|
|
bcoeffs = [my .- mX*bcoeffs; bcoeffs];
|
|
end
|
|
|
|
lambda_min, lambda_min_ind = findmin(rmsecv)
|
|
lambda_min_ind = lambda_min_ind[1]
|
|
b_lambda_min = bcoeffs[:,lambda_min_ind]
|
|
|
|
return b_lambda_min, rmsecv, lambda_min, lambda_min_ind
|
|
end
|
|
|
|
"""
|
|
function plegendre(d, p)
|
|
|
|
Calculates orthonormal Legendre polynomials using a QR factorisation.
|
|
Inputs:
|
|
- d : polynomial degree
|
|
- p : size of vector
|
|
|
|
Outputs:
|
|
- Q : (d+1) x p matrix with basis
|
|
- R : matrix from QR-factorisation
|
|
"""
|
|
function plegendre(d, p)
|
|
|
|
P = ones(p, d+1);
|
|
x = (-1:2/(p-1):1)';
|
|
for k in 1:d
|
|
P[:,k+1] = x.^k;
|
|
end
|
|
|
|
factorisation = qr(P);
|
|
Q = Matrix(factorisation.Q)';
|
|
R = Matrix(factorisation.R);
|
|
|
|
return Q, R
|
|
end
|