SpectralDP(Correlation)

 function [A,Corr] = SpectralDP(Correlation)

Corr_Mat=Correlation;
n=size(Corr_Mat,1);
Corr_Mat=Corr_Mat.*not(eye(n))+eye(n); 

%Spectral Decomposition
[Eig_Vec, ~]=eig(Corr_Mat); 
Eig_Val_Vector=eig(Corr_Mat); 
PosTest=Eig_Val_Vector >0; 
PosTest_n = Eig_Val_Vector <= 0;

% Use vectorized indexing
Eig_Val_Negative = Eig_Val_Vector;
Eig_Val_Positive = Eig_Val_Vector;
P_Ceiling_Idx= find(PosTest);
P_Floor_Idx = find(PosTest_n);
epsilon=0.0001;

% Break down negative and 0 and adjust - avoid for loop
Eig_Val_Negative(P_Floor_Idx) = Eig_Val_Negative(P_Floor_Idx)+ epsilon; 
Neg = Eig_Val_Negative(P_Floor_Idx).*(PosTest_n(P_Floor_Idx));

% Break down positive, and adjust to avoid smaller eigenvalues having
% a disproportionate influence.
Eig_Val_Positive(P_Ceiling_Idx) = Eig_Val_Positive(P_Ceiling_Idx)+epsilon;
Pos = Eig_Val_Positive(P_Ceiling_Idx).* (PosTest(P_Ceiling_Idx));

% Concatenate positive and non-negative
Eig_Val_Vector=[Neg; Pos];
Eig_Val_Mat=[Neg; Pos]*ones(1, n).*eye(n); 

% Re-construct new matrix and Cholesky
Compile = repmat(sum(Eig_Vec.*Eig_Vec.*repmat...
    (Eig_Val_Vector', n, 1), 2).^(-1), 1, n).*eye(n); 
L=(Compile.^0.5)*Eig_Vec*(Eig_Val_Mat.^0.5);
Corr=L*L';
A=chol(Corr);