docs: def asalsan(X, rank, **kwargs):
      """
      ASALSAN algorithm to compute the three-way DEDICOM decomposition
      of a tensor
  
      See
      ---
      .. [1] Brett W. Bader, Richard A. Harshman, Tamara G. Kolda
         "Temporal analysis of semantic graphs using ASALSAN"
         7th International Conference on Data Mining, 2007
  
      .. [2] Brett W. Bader, Richard A. Harshman, Tamara G. Kolda
         "Temporal analysis of Social Networks using Three-way DEDICOM"
         Technical Report, 2006
      """
      # init options
      ainit = kwargs.pop('init', _DEF_INIT)
      proj = kwargs.pop('proj', _DEF_PROJ)
      maxIter = kwargs.pop('maxIter', _DEF_MAXITER)
      conv = kwargs.pop('conv', _DEF_CONV)
      nne = kwargs.pop('nne', _DEF_NNE)
      optfunc = kwargs.pop('optfunc', _DEF_OPTFUNC)
      if not len(kwargs) == 0:
          raise BaseException('Unknown keywords (%s)' % (kwargs.keys()))
  
      # init starting points
      D = ones((len(X), rank))
      sz = X[0].shape
      n = sz[0]
      R = rand(rank, rank)
      if ainit == 'random':
          A = rand(n, rank)
      elif ainit == 'nvecs':
          S = zeros((n, n))
          T = zeros((n, n))
          for i in range(len(X)):
              T = X[i]
              S = S + T + T.T
          evals, A = eigsh(S, rank)
          if nne > 0:
              A[A < 0] = 0
          if proj:
              Q, A2 = qr(A)
              X2 = __projectSlices(X, Q)
              R = __updateR(X2, A2, D, R, nne)
          else:
              R = __updateR(X, A, D, R, nne)
      elif isinstance(ainit, np.ndarray):
          A = ainit
      else:
          raise 'Unknown init option ("%s")' % ainit
  
      # perform decomposition
      if issparse(X[0]):
          normX = [norm(M.data) ** 2 for M in X]
          Xflat = [M.tolil().reshape((1, prod(M.shape))).tocsr() for M in X]
      else:
          normX = [norm(M) ** 2 for M in X]
          Xflat = [M.flatten() for M in X]
      M = zeros((n, n))
      normXSum = sum(normX)
      #normX = norm(X)**2
      fit = fitold = f = fitchange = 0
      exectimes = []
      for iters in xrange(maxIter):
          tic = time.clock()
          fitold = fit
          A = __updateA(X, A, D, R, nne)
          if proj:
              Q, A2 = qr(A)
              X2 = __projectSlices(X, Q)
              R = __updateR(X2, A2, D, R, nne)
              D, f = __updateD(X2, A2, D, R, nne, optfunc)
          else:
              R = __updateR(X, A, D, R, nne)
              D, f = __updateD(X, A, D, R, nne, optfunc)
  
          # compute fit
          f = 0
          for i in xrange(len(X)):
              AD = dot(A, diag(D[i, :]))
              M = dot(dot(AD, R), AD.T)
              f += normX[i] + norm(M) ** 2 - 2 * Xflat[i].dot(M.flatten())
          f *= 0.5
          fit = 1 - (f / normXSum)
          fitchange = abs(fitold - fit)
  
          exectimes.append(time.clock() - tic)
  
          # print iter info when debugging is enabled
          _log.debug('[%3d] fit: %.5f | delta: %7.1e | secs: %.5f' % (
              iters, fit, fitchange, exectimes[-1]
          ))
  
          if iters > 1 and fitchange < conv:
              break
      return A, R, D, fit, iters, array(exectimes)