diff --git a/src/TR.jl b/src/TR.jl
index 21fb061..e5a0876 100644
--- a/src/TR.jl
+++ b/src/TR.jl
@@ -62,16 +62,28 @@ end
 
 function regularizationMatrix(p::Int64; regType="L2", regParam1=0, regParam2=1e-14)
     
-if regType == "bc"
+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"
+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"
+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)^(k-1) * (1-(regParam1+1)/k) * C[k-1];
+    end
+
+    regMat = zeros(p,p);
+
+    for i in 1:p
+        regMat[i:end, i] = regMat[1:end-i+1];
+    end
 end
 
 return regMat