code: def dkron(A,dA,B,dB, operation='prod'):
      
  
      if dA is None:
          dA_param_num = 0
          dA = np.zeros((A.shape[0], A.shape[1],1))
      else:
          dA_param_num = dA.shape[2]
  
      if dB is None:
          dB_param_num = 0
          dB = np.zeros((B.shape[0], B.shape[1],1))
      else:
          dB_param_num = dB.shape[2]
  
      # Space allocation for derivative matrix
      dC = np.zeros((A.shape[0]*B.shape[0], A.shape[1]*B.shape[1], dA_param_num +  dB_param_num))
  
      for k in range(dA_param_num):
          if operation == 'prod':
              dC[:,:,k] = np.kron(dA[:,:,k],B);
          else:
              dC[:,:,k] = np.kron(dA[:,:,k],np.eye( B.shape[0] ))
  
      for k in range(dB_param_num):
          if operation == 'prod':
              dC[:,:,dA_param_num+k] = np.kron(A,dB[:,:,k])
          else:
              dC[:,:,dA_param_num+k] = np.kron(np.eye( A.shape[0] ),dB[:,:,k])
  
      return dC
docs: Function computes the derivative of Kronecker product A*B
      (or Kronecker sum A+B).
  
      Input:
      -----------------------
  
      A: 2D matrix
          Some matrix
      dA: 3D (or 2D matrix)
          Derivarives of A
      B: 2D matrix
          Some matrix
      dB: 3D (or 2D matrix)
          Derivarives of B
  
      operation: str 'prod' or 'sum'
          Which operation is considered. If the operation is 'sum' it is assumed
          that A and are square matrices.s
  
      Output:
          dC: 3D matrix
          Derivative of Kronecker product A*B (or Kronecker sum A+B)