title: Python Implementation of K-SVD algorithm
body: <p>I am currently trying to apply a K-SVD algorithm using Python that is based on the following paper over here: <a href="http://www.cs.technion.ac.il/~elad/publications/journals/2004/32_KSVD_IEEE_TSP.pdf" rel="nofollow noreferrer">http://www.cs.technion.ac.il/~elad/publications/journals/2004/32_KSVD_IEEE_TSP.pdf</a> . So far I have implemented the code below and it seems to be working quite fine in terms of creating a good dictionary for sparse representations. The code that I created is based loosely on an implementation that exists in Python here, <a href="https://github.com/nel215/ksvd/blob/master/ksvd/__init__.py" rel="nofollow noreferrer">https://github.com/nel215/ksvd/blob/master/ksvd/<strong>init</strong>.py</a>. However I modified as I am using the K-svd to perform dictionary learning for a different kind of signals compared to the one normal K-svd is supposed to be used for (I am attempting to create a sparse dictionary for machine vibration signals instead of normal image signals). My code can be seen below:</p>
  
  <pre><code>def Phi_designer_k_svd(Phi_list, y_list, index):
      '''
      # y_list is expected to have size of k x y_freq (k rows of y_freq sized horizontal signals)
      # mp_process in here is a self-implemented OMP function, loosely following
      # an implementation here: https://github.com/davebiagioni/pyomp/blob/master/omp.py
      '''
      # (mxn) = size of Phi
      Phi = Phi_list[index]
  
      for i in range(0,1000):
          Phi_old = np.zeros(Phi.shape)
          Phi_old[:, :] = Phi[:, :]
          x_mp = np.zeros((n,k))
  
          #for every column of x_mp, calculate the sparse representation of the column
          #(find the representation of x_mp that would give minimum solution for ||y - Phi*x_mp||)
          for j in range(0, k): 
              # find approximation of x signal
              x_mp[:,j], _, _ = mp_process(Phi, y_cur[:,j], ncoef=sparsity, verbose=vbose)
  
          '''for every t-th atom in Phi...
          update the dictionary atoms so that it minimizes errors obtained from compressed x_mp
          '''
          for t in range(0, n):
  
              #Choose the COLUMN indexes in the t-th ROW of x_mp that is NON zero!
              #(synonymous with picking signals (which column) of x contributed 
              # directly by the t-th atom of the dictionary)
  
              I = x_mp[t] != 0
  
              #if there are no contributions made by this atom for x_mp, then continue to the next atom.
              if np.sum(I) == 0:
                  continue
  
  
              '''
              # only columns containing nonzero elements from t-th row of x_mp is used (showing indices using t-th atom in dict), 
              # rest are ignored
              '''
              x_copy = x_mp[:, I]
  
              #zero the t-th row as it will not be used in actual calculation
              x_copy[t] = 0
  
              #create a copy of Phi with the value of the t-th atom zeroed (to disable contribution from the t-th atom of Phi)
              copy = np.zeros(Phi.shape)
              copy[:] = Phi[:]
              copy[:,t] = 0
  
              #calculate error produced from contribution of t-th atom only (thus ignoring the rest of the zero elements in initial x_mp.
              error = y_cur[:,I] - np.dot(copy, x_copy)
  
  
              #produce a SVD decomp of the obtained error matrix
              U,s,V = np.linalg.svd(error, full_matrices=True)
  
              Phi[:, t] = U[:, 0]
  
              '''
              #update only the picked non-zero elements of x_mp (as previously mentioned) to be updated. 
              #(sizes of s and V should have already matched this indices group as well)
              '''
              x_mp[t, I] = s[0] * V[:,0]
  
  
          previous_norm = l2_norm(Phi_old)
          detected_norm = l2_norm(Phi)
  
          norm_diff = previous_norm - detected_norm
  
          #**** Convergence condition. Not sure correct or not based on paper.
          if abs(norm_diff) &lt; tol:
              break
  
      Phi_list[index] = Phi
      return
  </code></pre>
  
  <p>However I found that this code is only able to create a dictionary that works well for data training set only, and the compression does not work for any new sets of data. I am suspicious that I may have not implemented the stopping criteria from the paper correctly (the one that I marked with '****' in the code above. Last line in the for loop), as I am not too sure my self what would the appropriate stopping criteria for the algorithm be. Has anybody ever implemented their own version of the K-SVD in python before? If so I would really appreciate if you can help me check the validity of my implemented stopping criteria, and perhaps any kind of extra tips that I can use to ensure better match between the created dictionary and any new set of signals.</p>