implement SVD for python examples

python/pseudoinverse.py

@@ -0,0 +1,33 @@
+import numpy as np
+import scipy as sp
+
+"""
+This algorithm solves the Thikonov regularized least squares problem 
+using the singular value decomposition of A.
+
+# Arguments
+
+* `U, Sigm, V`: Singular value decomposition of A
+* `b::Vector`: Measurement vector b
+* `lambd::Float64`: The regularization parameter, relative to the matrix trace
+* `enforceReal::Bool`: Enable projection of solution on real plane during iteration
+* `enforcePositive::Bool`: Enable projection of solution onto positive halfplane during iteration
+"""
+def pseudoinverse(U, Sigm, V, b, lambd, enforceReal, enforcePositive):
+  # perform regularization
+  D = np.zeros(np.size(Sigm))
+  for i in range(np.size(Sigm)):
+    sigmi = Sigm[i]
+    D[i] = sigmi/(sigmi**2+lambd**2)
+
+  # calculate pseudoinverse
+  tmp = np.dot(U.conjugate().transpose(),b[:]) #  conjugate transpose
+  tmp =  tmp*D
+  c = np.dot(V.conjugate().transpose(),tmp) # not transposed
+
+  if enforceReal and np.iscomplexobj(c):
+    c.imag = 0
+  if enforcePositive:
+    c = c * (c.real > 0)
+
+  return c

python/reco.py

@@ -1,4 +1,5 @@
 from kaczmarzReg import *
+from pseudoinverse import *
 from pylab import *
 import h5py
 import urllib
@@ -45,10 +46,24 @@
 # reconstruct
 c = kaczmarzReg(S,u,1,1e6,False,True,True)
 
+
+# reconstruct using signular value decomposition
+U, Sigm, V = svd(S, full_matrices=False)
+csvd = pseudoinverse(U, Sigm, V, u, 5e3, True, True)
+
 # reshape into an image
 N = fSM['/calibration/size'][:]
 c = reshape(c,(N[0],N[1]))
+csvd = reshape(csvd,(N[0],N[1]))
 
+# plot kaczmarz reconstruction
+figure()
 gray()
-imshow(real(c))
+imshow(real(transpose(c)), interpolation="None")
+
+# plot pseudoinverse reconstruction
+figure()
+gray()
+imshow(real(transpose(csvd)), interpolation="None")
+
 show()