diff --git a/computeFlattening.m b/computeFlattening.m
index b8f2221..bb7e2ff 100644
--- a/computeFlattening.m
+++ b/computeFlattening.m
@@ -11,7 +11,15 @@
 M=[L A'; A sparse(n_eq,n_eq)];
 rhs=[zeros(n_vars,1); b];
 warning('off','MATLAB:nearlySingularMatrix')
-x_lambda = M \ rhs;
+
+% Use Tikhonov Regularization for calculating x_lambda
+% Reference: Tikhonov, A. N., & Arsenin, V. Y. (1977). Solutions of Ill-Posed Problems. Winston & Sons.
+% addresses rank deficiency in the original M matrix
+% specify alpha: 
+% - small enough to meet error requirements
+% - large enough to make eignvalues not zero 
+alpha = 1e-10;
+x_lambda = (M' * M + alpha * speye(size(M,2))) \ (M' * rhs);
 warning('on','MATLAB:nearlySingularMatrix')
 e=max(abs(M*x_lambda-rhs));
 fprintf('error: %e\n',e);