docs: def low_rank_approximation(A, r):
      """ Returns approximation of A of rank r in least-squares sense."""
      try:
          u, s, v = np.linalg.svd(A, full_matrices=False)
      except np.linalg.LinAlgError as e:
          print('Matrix:', A)
          print('Matrix rank:', np.linalg.matrix_rank(A))
          raise
  
      Ar = np.zeros((len(u), len(v)), dtype=u.dtype)
      buf_ = np.empty_like(Ar)
      sc_vec_ = np.empty((v.shape[1],), dtype=v.dtype)
      for i in range(r):
          np.multiply(v[i], s[i], out=sc_vec_)
          np.outer(u[:, i], sc_vec_, out=buf_)
          Ar += buf_
      return Ar