diff --git a/factorize.m b/factorize.m
new file mode 100644
index 0000000000000000000000000000000000000000..ee801cddee41877b29848445543b999c9e5b4b9a
--- /dev/null
+++ b/factorize.m
@@ -0,0 +1,14 @@
+function [B,C,U,S,V] = factorize(M,k)
+%FACTORIZE Summary of this function goes here
+%   Detailed explanation goes here
+
+Mask = isnan(M);
+M(Mask) = 0;
+
+[U,S,V] = utils.svd_nsq(M);
+B = U(:,1:k)*S(1:k,1:k);
+C = V(:,1:k)';
+
+B(any(Mask,2),:) = nan;
+
+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