code: def sqize_kernel_low_rank_fast(
      K,
      L,
      theta=1.0,
      normalize=False,
      rankbound=float("inf"),
  ):
      """Compute the sequential kernel from kernel matrix, with low-rank approximation.
  
      Vectorized across series, and using (faster) low-rank approximation across series.
  
      Parameters
      ----------
      K : 3D np.ndarray
          containing joint low-rank factors (symmetric singular factor)
              1st index counts sequences
              2nd index counts time
              3rd index counts features
          so K[m,:,:] is the mth factor,
          and K[m,:,:] x K[m,:,:]^t is the kernel matrix of the mth factor
      L : an integer >= 1, representing the level of truncation
      theta : float, optional, default=1.0
          a positive scaling factor for the levels, i-th level is scaled by theta^i
      normalize : bool, optional, default = False
          whether the output kernel matrix is normalized
          if True, sums and cumsums are divided by prod(K.shape)
      rankbound : int, optional, default = infinity
          a hard threshold for the rank of the level matrices
  
      Returns
      -------
      np.ndarray of shape (m, r), where r = min(rankbound, K.shape[2])
          a matrix R such that R*R^t is the sequential kernel matrix
          R*R^t[i,j] is sequential kernel (low-rank) between i-th and j-th sequence in K
      """
      if normalize:
          Ksize = K.shape[0]
          B = np.ones([Ksize, 1, 1])
          R = np.ones([Ksize, 1])
  
          for _ in range(L):
              P = np.sqrt(theta) * hadamard_low_rank_batch(K, B) / Ksize
              B = cumsum_shift_mult(P, [1])
  
              if rankbound < B.shape[2]:
                  # B = rankreduce_batch(B, rankbound)
                  permut = np.sort(
                      np.random.permutation(range(B.shape[2]))[
                          range(rankbound)
                      ]
                  )
                  B = B[:, :, permut]
  
              R = np.concatenate((R, np.sum(B, axis=1)), axis=1) / (
                  np.sqrt(1 + theta)
              )
  
          return R
  
      else:
          Ksize = K.shape[0]
          B = np.ones([Ksize, 1, 1])
          R = np.ones([Ksize, 1])
  
          for _ in range(L):
              # todo: execute only if rank is lower than rankbound
              #       reduce to rank
              P = np.sqrt(theta) * hadamard_low_rank_batch(K, B)
              B = cumsum_shift_mult(P, [1])
  
              if rankbound < B.shape[2]:
                  # B = rankreduce_batch(B, rankbound)
                  permut = np.sort(
                      np.random.permutation(range(B.shape[2]))[
                          range(rankbound)
                      ]
                  )
                  B = B[:, :, permut]
  
              R = np.concatenate((R, np.sum(B, axis=1)), axis=1)
  
          return R
  
imports: import collections
  from functools import partial
  import numpy
  from scipy.sparse.linalg import svds
  from sklearn.base import BaseEstimator, TransformerMixin
  from sktime.dists_kernels.base import BasePairwiseTransformerPanel
label_desc: Use image node_id 29 for calling a global function with example usage: sqize_kernel_low_rank_fast(K, L, theta, normalize, rankbound) and returns: R, R