diff --git a/factorize.m b/factorize.m
new file mode 100644
index 0000000000000000000000000000000000000000..62a64d18602dd77069e9016ee699e257e8fe8e2c
--- /dev/null
+++ b/factorize.m
@@ -0,0 +1,27 @@
+function [LFB,C,U,S,V] = factorize(LFRef,varargin)
+%FACTORIZE Summary of this function goes here
+%   Detailed explanation goes here
+
+sz = size(LFRef);
+Mask = isnan(LFRef); LFRef(Mask) = 0;
+vecResh = @(x)   reshape(x,size(x,1)*size(x,2),[])';
+boxResh = @(x,sz) reshape(x',sz);
+
+LFRefvec = vecResh(LFRef);
+Maskvec = vecResh(Mask);
+
+kList = 1:size(LFRefvec,2);
+if nargin>1
+    kList = varargin{1};
+end
+
+[U,S,V] = svd_nsq(LFRefvec);
+B = U*S;
+C = V';
+
+B(Maskvec) = NaN;
+B = B(:,kList);
+C = C(kList,:);
+
+LFB = boxResh(B,[numel(kList),1,sz(3:end)]);
+end
\ No newline at end of file
diff --git a/svd_nsq.m b/svd_nsq.m
new file mode 100644
index 0000000000000000000000000000000000000000..f7521ce7902916f81041bf8aaf2ec559ad06264a
--- /dev/null
+++ b/svd_nsq.m
@@ -0,0 +1,28 @@
+function [U,S,V] = svd_nsq(Z)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%               SVD accelerated for non-square matrices
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+[m, n] = size(Z);
+
+if 2*m < n
+    ZZT = Z*Z';
+    [U, S2, Ua] = svd(ZZT);
+    S = sqrt(S2);
+    Sinv = diag(1./diag(S));
+    V = Z' * U * Sinv;
+    return;
+end
+if m > 2*n
+    ZTZ = Z'*Z;
+    [Va, S2, V] = svd(ZTZ);
+    S = sqrt(S2);
+    Sinv = diag(1./diag(S));
+    U = Z * V * Sinv;
+    return;
+end
+
+if(m==n)
+    [U,S,V] = svd(Z);
+else
+    [U,S,V] = svd(Z,'econ');
+end
\ No newline at end of file