Skip to content

API

Table of Contents

Optimal weights summing to 1 and WCE

opt_weights_sum_1 

opt_weights_sum_1(x, kernel, dims=2, eps=1e-12, return_wce2_p_kmean=True, use_pcg=False, pcg_kwargs={'verbose': False})

Compute the optimal weights for a kernel quadrature rule.

Parameters:

Name Type Description Default
x Tensor

A batch of samples of shape (..., N, *m) where len(m)==dims.

 required
kernel callable

Kernel function.

 required
dims int

See description of x.

2
eps float

Regularization parameter for numerical stability.

 1e-12
return_wce2_p_kmean bool

If True, return \(\mathrm{WCE}^2+\gamma = \boldsymbol{w}^\intercal \mathsf{K} \boldsymbol{w}\) where \(\gamma = \int_\mathcal{M} K(\boldsymbol{x},\boldsymbol{y}) \mu(\mathrm{d} \boldsymbol{x})\) is assumed to be constant as a function of \(\boldsymbol{x}\).

True
use_pcg bool

If True, use preconditioned conjugate gradient (PCG) to solve for weights. If use_pcg an int is passed in, then PCG is used whenever N>=use_pcg.

 False
pcg_kwargs dict

keyword arguments to pass to agsutil.pcg.

 {'verbose': False}

Returns:

Name Type Description
w Tensor

Optimal weights of shape (..., N).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> kernel = KernelMatern(nu=1/2)
>>> x = rand_flag_R(3,7,seed=7)
>>> x.shape
torch.Size([3, 7, 7])
>>> w,wce2_p_kmean = opt_weights_sum_1(x,kernel)
>>> w.shape 
torch.Size([3])
>>> w.sum() 
tensor(1.)
>>> w
tensor([0.3506, 0.3218, 0.3276])
>>> wce2_p_kmean
tensor(0.7027)

>>> x2 = rand_flag_R(3,7,seed=11)
>>> xfull = torch.stack([x,x2],dim=0)
>>> xfull.shape 
torch.Size([2, 3, 7, 7])
>>> w_full,wce2_p_kmean_full = opt_weights_sum_1(xfull,kernel)
>>> w_full.shape 
torch.Size([2, 3])
>>> w_full.sum(-1) 
tensor([1., 1.])
>>> w_full
tensor([[0.3506, 0.3218, 0.3276],
        [0.3015, 0.3055, 0.3931]])
>>> wce2_p_kmean_full.shape 
torch.Size([2])
>>> wce2_p_kmean_full
tensor([0.7027, 0.7009])

>>> w_full_pcg,wce2_p_kmean_full_pcg = opt_weights_sum_1(xfull,kernel,use_pcg=True)
>>> torch.allclose(w_full,w_full_pcg)
True
>>> torch.allclose(wce2_p_kmean_full,wce2_p_kmean_full_pcg)
True
Source code in rmatqmcpy/opt.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def opt_weights_sum_1(x, kernel, dims=2, eps=1e-12, return_wce2_p_kmean=True, use_pcg=False, pcg_kwargs={"verbose":False}):
    r"""
    Compute the optimal weights for a kernel quadrature rule.

    Args:
        x (torch.Tensor): A batch of samples of shape `(..., N, *m)` where `len(m)==dims`.
        kernel (callable): Kernel function.
        dims (int): See description of `x`. 
        eps (float, optional): Regularization parameter for numerical stability.
        return_wce2_p_kmean (bool): If True, return 
            $\mathrm{WCE}^2+\gamma = \boldsymbol{w}^\intercal \mathsf{K} \boldsymbol{w}$ 
            where $\gamma = \int_\mathcal{M} K(\boldsymbol{x},\boldsymbol{y}) \mu(\mathrm{d} \boldsymbol{x})$ 
            is assumed to be constant as a function of $\boldsymbol{x}$. 
        use_pcg (bool): If True, use preconditioned conjugate gradient (PCG) to solve for weights. 
            If `use_pcg` an int is passed in, then PCG is used whenever `N>=use_pcg`.
        pcg_kwargs (dict): keyword arguments to pass to `agsutil.pcg`.

    Returns:
        w (torch.Tensor): Optimal weights of shape `(..., N)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> kernel = KernelMatern(nu=1/2)
        >>> x = rand_flag_R(3,7,seed=7)
        >>> x.shape
        torch.Size([3, 7, 7])
        >>> w,wce2_p_kmean = opt_weights_sum_1(x,kernel)
        >>> w.shape 
        torch.Size([3])
        >>> w.sum() 
        tensor(1.)
        >>> w
        tensor([0.3506, 0.3218, 0.3276])
        >>> wce2_p_kmean
        tensor(0.7027)

        >>> x2 = rand_flag_R(3,7,seed=11)
        >>> xfull = torch.stack([x,x2],dim=0)
        >>> xfull.shape 
        torch.Size([2, 3, 7, 7])
        >>> w_full,wce2_p_kmean_full = opt_weights_sum_1(xfull,kernel)
        >>> w_full.shape 
        torch.Size([2, 3])
        >>> w_full.sum(-1) 
        tensor([1., 1.])
        >>> w_full
        tensor([[0.3506, 0.3218, 0.3276],
                [0.3015, 0.3055, 0.3931]])
        >>> wce2_p_kmean_full.shape 
        torch.Size([2])
        >>> wce2_p_kmean_full
        tensor([0.7027, 0.7009])

        >>> w_full_pcg,wce2_p_kmean_full_pcg = opt_weights_sum_1(xfull,kernel,use_pcg=True)
        >>> torch.allclose(w_full,w_full_pcg)
        True
        >>> torch.allclose(wce2_p_kmean_full,wce2_p_kmean_full_pcg)
        True
    """
    assert x.ndim>=(dims+1)
    kmat = kernel(x.unsqueeze(-dims-1),x.unsqueeze(-dims-2),dim=tuple(j for j in range(-dims,0,1))) # (...,N,N)
    N = kmat.size(-1)
    b = torch.ones(kmat.shape[:-1],device=x.device) # (...,N)
    if not isinstance(use_pcg,bool):
        assert isinstance(use_pcg,int) 
        use_pcg = N>=use_pcg
    assert isinstance(use_pcg,bool)
    kmat_noisy = kmat+eps*torch.eye(N,device=x.device)
    if use_pcg:
        w = agsutil.pcg(
            kmat_noisy,
            b[...,None],
            **pcg_kwargs,
        )[...,0]
    else:
        L = torch.linalg.cholesky(kmat+eps*torch.eye(N,device=x.device),upper=False) # (...,N,N)
        w = torch.cholesky_solve(b[...,None],L,upper=False)[...,0] # (...,N)
    w = w/w.sum(-1,keepdim=True)
    if return_wce2_p_kmean:
        wce_squared_plus_kernel_integral = torch.einsum("...i,...ij,...j->...",w,kmat,w) # (...,)
        return w,wce_squared_plus_kernel_integral
    else:
        return w

Matern and polynomial kernels

KernelMatern 

KernelMatern(nu=5 / 2, sigma2=1, rho=1)

Bases: object

Matern chordal kernel.

Parameters:

Name Type Description Default
nu float

Smoothness parameter. Must be in [1/2, 3/2, 5/2, 7/2, 9/2, 11/2, 13/2, np.inf].

 5 / 2
sigma2 float

Variance parameter.

 1
rho float

Length-scale parameter.

 1

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> kernel = KernelMatern(nu=1/2)

>>> x = rand_flag_R(3,7,seed=7)
>>> x.shape 
torch.Size([3, 7, 7])
>>> k = kernel(x[:,None,:,:],x[None,:,:,:])
>>> k.shape 
torch.Size([3, 3])
>>> k
tensor([[1.0000, 0.5458, 0.5388],
        [0.5458, 1.0000, 0.5787],
        [0.5388, 0.5787, 1.0000]])

>>> for rand_manifold in [
...         rand_flag_R,
...         rand_flag_C,
...         rand_flag_H,
...         rand_Stiefel_R,
...         rand_Stiefel_C,
...         rand_Stiefel_H,
...         rand_Gr_R,
...         rand_Gr_C,
...         rand_Gr_H,
...         rand_LGr_R,
...         rand_LGr_C,
...         rand_LGr_H,]:
...     x = rand_manifold(3,7,seed=7)
...     for nu in [1/2,3/2,5/2,7/2,9/2,11/2,13/2,np.inf]:
...         kernel = KernelMatern(nu=nu)
...         k = kernel(x[:,None,:,:],x[None,:,:,:])
...         assert k.shape==(3,3)
...         assert not k.isnan().any()
...         assert k.isfinite().all()
...         assert not k.is_complex()
Source code in rmatqmcpy/kernels.py 
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
def __init__(self, nu=5/2, sigma2=1, rho=1):
    self.nu = nu
    self.sigma2 = sigma2 
    self.rho = rho 
    assert self.nu in self.SUPPORTED_NU, "nu should be in %s"%str(self.SUPPORTED_NU)
    if self.nu==1/2: 
        self.c = torch.tensor([1])
        self.nu_str = r"$\nu = 1/2$"
    elif self.nu==3/2:
        self.c = torch.tensor([1,np.sqrt(3)])
        self.nu_str = r"$\nu = 3/2$"
    elif self.nu==5/2:
        self.c = torch.tensor([1,np.sqrt(5),5/3])
        self.nu_str = r"$\nu = 5/2$"
    elif self.nu==7/2:
        self.c = torch.tensor([1,np.sqrt(7),14/5,7*np.sqrt(7)/15])
        self.nu_str = r"$\nu = 7/2$"
    elif self.nu==9/2:
        self.c = torch.tensor([1,3,27/7,18/7,27/35])
        self.nu_str = r"$\nu = 9/2$"
    elif self.nu==11/2:
        self.c = torch.tensor([1,np.sqrt(11),44/9,11*np.sqrt(11)/9,121/63,121*np.sqrt(11)/945])
        self.nu_str = r"$\nu = 11/2$"
    elif self.nu==13/2:
        self.c = torch.tensor([1,np.sqrt(13),65/11,52*np.sqrt(13)/33,338/99,169*np.sqrt(13)/495,2197/10395])
        self.nu_str = r"$\nu = 13/2$"
    elif self.nu==np.inf:
        self.nu_str = r"$\nu = \infty$"

__call__ 

__call__(x, y, dim=(-2, -1))

Compute the Matern kernel between two batches of tensors.

Parameters:

Name Type Description Default
x Tensor

First batch of tensors.

 required
y Tensor

Second batch of tensors.

 required
dim tuple

Dimensions over which to take the Frobenius distance..

 (-2, -1)

Returns:

Name Type Description
k Tensor

The computed Matern kernel values.
Source code in rmatqmcpy/kernels.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def __call__(self, x, y, dim=(-2,-1)):
    r"""
    Compute the Matern kernel between two batches of tensors.

    Args:
        x (torch.Tensor): First batch of tensors.
        y (torch.Tensor): Second batch of tensors.
        dim (tuple, optional): Dimensions over which to take the Frobenius distance..

    Returns:
        k (torch.Tensor): The computed Matern kernel values.
    """
    frob_dist = torch.linalg.norm(x-y,dim=dim).real
    rscaled = frob_dist/self.rho
    if self.nu==np.inf:
        k = torch.exp(-(rscaled)**2/2)
    else:
        exp_term = torch.exp(-np.sqrt(2*self.nu)*rscaled) 
        powers = torch.arange(len(self.c),dtype=x.real.dtype,device=x.device)
        poly_term = (self.c.to(x.device)*rscaled[...,None]**powers).sum(-1)
        k = exp_term*poly_term
    return self.sigma2*k

KernelPoly 

KernelPoly(t=2, c=1)

Bases: object

Polynomial kernel.

Args: r (int): Degree of the polynomial. c (float): Constant offset.

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> kernel = KernelPoly()

>>> x = rand_flag_R(3,7,seed=7)
>>> x.shape 
torch.Size([3, 7, 7])
>>> k = kernel(x[:,None,:,:],x[None,:,:,:])
>>> k.shape 
torch.Size([3, 3])
>>> k
tensor([[4.0000, 3.3004, 3.2718],
        [3.3004, 4.0000, 3.4239],
        [3.2718, 3.4239, 4.0000]])

>>> for rand_manifold in [
...         rand_flag_R,
...         rand_flag_C,
...         rand_flag_H,
...         rand_Stiefel_R,
...         rand_Stiefel_C,
...         rand_Stiefel_H,
...         rand_Gr_R,
...         rand_Gr_C,
...         rand_Gr_H,
...         rand_LGr_R,
...         rand_LGr_C,
...         rand_LGr_H,]:
...     x = rand_manifold(3,7,seed=7)
...     k = kernel(x[:,None,:,:],x[None,:,:,:])
...     assert k.shape==(3,3)
...     assert not k.isnan().any()
...     assert k.isfinite().all()
...     assert not k.is_complex()
Source code in rmatqmcpy/kernels.py 
163
164
165
166
167
168
def __init__(self, t=2, c=1):
    assert t>=1 
    assert t%1==0
    assert c>0
    self.t = t
    self.c = c

__call__ 

__call__(x, y, dim=(-2, -1))

Compute the polynomial kernel between two batches of tensors.

Parameters:

Name Type Description Default
x Tensor

First batch of tensors.

 required
y Tensor

Second batch of tensors.

 required
dim tuple

Dimensions over which to take the Frobenius inner product.

 (-2, -1)

Returns:

Name Type Description
k Tensor

The computed polynomial kernel values.
Source code in rmatqmcpy/kernels.py 
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
def __call__(self, x, y, dim=(-2,-1)):
    r"""
    Compute the polynomial kernel between two batches of tensors.

    Args:
        x (torch.Tensor): First batch of tensors.
        y (torch.Tensor): Second batch of tensors.
        dim (tuple, optional): Dimensions over which to take the Frobenius inner product.

    Returns:
        k (torch.Tensor): The computed polynomial kernel values.
    """
    frob_inner_prod = (x*y.conj()).sum(dim=dim).real
    v = (self.c+frob_inner_prod)**self.t
    return v

Random Stiefel, Grassmann, flag, and Lagrangian Grassmann matrices

rand_flag_R 

rand_flag_R(N, n, delta=None, seed=None, rand_On=None, qp_unif_gen=None, device=None)

Generate a batch of N random real flag matrices of size (N, n, n).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Dimension of the square matrix.

 required
delta Tensor

Eigenvalues of the flag matrices, defaults to torch.arange(1,n+1)/torch.linalg.norm(torch.arange(1,n+1)).

 None
seed int

Random seed for reproducibility.

 None
rand_On callable

Function to generate random orthogonal matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random real flag matrices of size (N, n, n).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_flag_R(2,3,seed=7)
>>> x.shape 
torch.Size([2, 3, 3])
>>> x
tensor([[[ 0.3107,  0.1198,  0.0819],
         [ 0.1198,  0.7127,  0.0915],
         [ 0.0819,  0.0915,  0.5802]],

        [[ 0.4100, -0.2117, -0.0450],
         [-0.2117,  0.5880,  0.1135],
         [-0.0450,  0.1135,  0.6057]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def rand_flag_R(N, n, delta=None, seed=None, rand_On=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random real flag matrices of size `(N, n, n)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Dimension of the square matrix.
        delta (torch.Tensor): Eigenvalues of the flag matrices, defaults to `torch.arange(1,n+1)/torch.linalg.norm(torch.arange(1,n+1))`.
        seed (int): Random seed for reproducibility.
        rand_On (callable): Function to generate random orthogonal matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random real flag matrices of size `(N, n, n)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_flag_R(2,3,seed=7)
        >>> x.shape 
        torch.Size([2, 3, 3])
        >>> x
        tensor([[[ 0.3107,  0.1198,  0.0819],
                 [ 0.1198,  0.7127,  0.0915],
                 [ 0.0819,  0.0915,  0.5802]],
        <BLANKLINE>
                [[ 0.4100, -0.2117, -0.0450],
                 [-0.2117,  0.5880,  0.1135],
                 [-0.0450,  0.1135,  0.6057]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if delta is None: 
        delta = torch.arange(1,n+1,device=device,dtype=torch.get_default_dtype())
        delta = delta/torch.linalg.norm(delta)
    assert delta.shape==(n,)
    assert not delta.is_complex()
    if rand_On is None: rand_On = rand_On_QR
    assert rand_On in [rand_On_QR,rand_On_eig]
    assert not torch.is_complex(delta)
    n = delta.size(-1)
    q = rand_On(N=N,n=n,k=n,seed=seed,device=delta.device,qp_unif_gen=qp_unif_gen)
    x = torch.einsum("...ij,...j,...kj->...ik",q,delta.to(q.dtype),q)
    return x

rand_flag_C 

rand_flag_C(N, n, delta=None, seed=None, rand_Un=None, qp_unif_gen=None, device=None)

Generate a batch of N random complex flag matrices of size (N, n, n).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Dimension of the square matrix.

 required
delta Tensor

Eigenvalues of the flag matrices, defaults to torch.arange(1,n+1)/torch.linalg.norm(torch.arange(1,n+1)).

 None
seed int

Random seed for reproducibility.

 None
rand_Un callable

Function to generate random unitary matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random complex flag matrices of size (N, n, n).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_flag_C(2,3,seed=7)
>>> x.shape 
torch.Size([2, 3, 3])
>>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
tensor([[[ 0.5330+0.0000j,  0.0552+0.1661j, -0.1884-0.0685j],
         [ 0.0552-0.1661j,  0.5332+0.0000j,  0.0186+0.0154j],
         [-0.1884+0.0685j,  0.0186-0.0154j,  0.5374+0.0000j]],

        [[ 0.4512+0.0000j, -0.1101+0.0152j,  0.1480-0.1704j],
         [-0.1101-0.0152j,  0.5788+0.0000j, -0.0483+0.0246j],
         [ 0.1480+0.1704j, -0.0483-0.0246j,  0.5736+0.0000j]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
def rand_flag_C(N, n, delta=None, seed=None, rand_Un=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random complex flag matrices of size `(N, n, n)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Dimension of the square matrix.
        delta (torch.Tensor): Eigenvalues of the flag matrices, defaults to `torch.arange(1,n+1)/torch.linalg.norm(torch.arange(1,n+1))`.
        seed (int): Random seed for reproducibility.
        rand_Un (callable): Function to generate random unitary matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random complex flag matrices of size `(N, n, n)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_flag_C(2,3,seed=7)
        >>> x.shape 
        torch.Size([2, 3, 3])
        >>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
        tensor([[[ 0.5330+0.0000j,  0.0552+0.1661j, -0.1884-0.0685j],
                 [ 0.0552-0.1661j,  0.5332+0.0000j,  0.0186+0.0154j],
                 [-0.1884+0.0685j,  0.0186-0.0154j,  0.5374+0.0000j]],
        <BLANKLINE>
                [[ 0.4512+0.0000j, -0.1101+0.0152j,  0.1480-0.1704j],
                 [-0.1101-0.0152j,  0.5788+0.0000j, -0.0483+0.0246j],
                 [ 0.1480+0.1704j, -0.0483-0.0246j,  0.5736+0.0000j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if delta is None: 
        delta = torch.arange(1,n+1,device=device,dtype=torch.get_default_dtype())
        delta = delta/torch.linalg.norm(delta)
    assert delta.shape==(n,)
    assert not delta.is_complex()
    if rand_Un is None: rand_Un = rand_Un_QR
    assert rand_Un in [rand_Un_QR,rand_Un_eig]
    assert not torch.is_complex(delta)
    n = delta.size(-1)
    q = rand_Un(N=N,n=n,k=n,seed=seed,device=delta.device,qp_unif_gen=qp_unif_gen)
    x = torch.einsum("...ij,...j,...kj->...ik",q,delta.to(q.dtype),q.conj())
    return x

rand_flag_H 

rand_flag_H(N, n, delta=None, seed=None, rand_Spn=None, qp_unif_gen=None, device=None)

Generate a batch of N random quaternionic flag matrices of size (N, 2n, 2n).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Dimension of the square matrix.

 required
delta Tensor

Eigenvalues of the flag matrices, defaults to torch.arange(1,n+1)/torch.linalg.norm(torch.arange(1,n+1)).

 None
seed int

Random seed for reproducibility.

 None
rand_Spn callable

Function to generate random symplectic unitary matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random quaternionic flag matrices of size (N, 2n, 2n).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_flag_H(2,3,seed=7)
>>> x.shape 
torch.Size([2, 6, 6])
>>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
tensor([[[ 0.5221+0.0000e+00j,  0.0063-1.2520e-01j,  0.0377-5.4600e-02j,
           0.0000+0.0000e+00j,  0.0484+4.0300e-02j,  0.0241+1.6270e-01j],
         [ 0.0063+1.2520e-01j,  0.5409+0.0000e+00j, -0.0698-1.0140e-01j,
          -0.0484-4.0300e-02j,  0.0000+0.0000e+00j, -0.0412+5.7600e-02j],
         [ 0.0377+5.4600e-02j, -0.0698+1.0140e-01j,  0.5406+0.0000e+00j,
          -0.0241-1.6270e-01j,  0.0412-5.7600e-02j,  0.0000+0.0000e+00j],
         [ 0.0000+0.0000e+00j, -0.0484+4.0300e-02j, -0.0241+1.6270e-01j,
           0.5221+0.0000e+00j,  0.0063+1.2520e-01j,  0.0377+5.4600e-02j],
         [ 0.0484-4.0300e-02j,  0.0000+0.0000e+00j,  0.0412+5.7600e-02j,
           0.0063-1.2520e-01j,  0.5409+0.0000e+00j, -0.0698+1.0140e-01j],
         [ 0.0241-1.6270e-01j, -0.0412-5.7600e-02j,  0.0000+0.0000e+00j,
           0.0377-5.4600e-02j, -0.0698-1.0140e-01j,  0.5406+0.0000e+00j]],

        [[ 0.4803+0.0000e+00j, -0.0923-1.2690e-01j,  0.1221+3.1900e-02j,
           0.0000+0.0000e+00j, -0.0532-1.0000e-04j,  0.0830-8.8000e-03j],
         [-0.0923+1.2690e-01j,  0.4942+0.0000e+00j,  0.0287+3.6500e-02j,
           0.0532+1.0000e-04j,  0.0000+0.0000e+00j, -0.0995+4.8100e-02j],
         [ 0.1221-3.1900e-02j,  0.0287-3.6500e-02j,  0.6290+0.0000e+00j,
          -0.0830+8.8000e-03j,  0.0995-4.8100e-02j,  0.0000+0.0000e+00j],
         [ 0.0000+0.0000e+00j,  0.0532-1.0000e-04j, -0.0830-8.8000e-03j,
           0.4803+0.0000e+00j, -0.0923+1.2690e-01j,  0.1221-3.1900e-02j],
         [-0.0532+1.0000e-04j,  0.0000+0.0000e+00j,  0.0995+4.8100e-02j,
          -0.0923-1.2690e-01j,  0.4942+0.0000e+00j,  0.0287-3.6500e-02j],
         [ 0.0830+8.8000e-03j, -0.0995-4.8100e-02j,  0.0000+0.0000e+00j,
           0.1221+3.1900e-02j,  0.0287+3.6500e-02j,  0.6290+0.0000e+00j]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
def rand_flag_H(N, n, delta=None, seed=None, rand_Spn=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random quaternionic flag matrices of size `(N, 2n, 2n)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Dimension of the square matrix.
        delta (torch.Tensor): Eigenvalues of the flag matrices, defaults to `torch.arange(1,n+1)/torch.linalg.norm(torch.arange(1,n+1))`.
        seed (int): Random seed for reproducibility.
        rand_Spn (callable): Function to generate random symplectic unitary matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random quaternionic flag matrices of size `(N, 2n, 2n)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_flag_H(2,3,seed=7)
        >>> x.shape 
        torch.Size([2, 6, 6])
        >>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
        tensor([[[ 0.5221+0.0000e+00j,  0.0063-1.2520e-01j,  0.0377-5.4600e-02j,
                   0.0000+0.0000e+00j,  0.0484+4.0300e-02j,  0.0241+1.6270e-01j],
                 [ 0.0063+1.2520e-01j,  0.5409+0.0000e+00j, -0.0698-1.0140e-01j,
                  -0.0484-4.0300e-02j,  0.0000+0.0000e+00j, -0.0412+5.7600e-02j],
                 [ 0.0377+5.4600e-02j, -0.0698+1.0140e-01j,  0.5406+0.0000e+00j,
                  -0.0241-1.6270e-01j,  0.0412-5.7600e-02j,  0.0000+0.0000e+00j],
                 [ 0.0000+0.0000e+00j, -0.0484+4.0300e-02j, -0.0241+1.6270e-01j,
                   0.5221+0.0000e+00j,  0.0063+1.2520e-01j,  0.0377+5.4600e-02j],
                 [ 0.0484-4.0300e-02j,  0.0000+0.0000e+00j,  0.0412+5.7600e-02j,
                   0.0063-1.2520e-01j,  0.5409+0.0000e+00j, -0.0698+1.0140e-01j],
                 [ 0.0241-1.6270e-01j, -0.0412-5.7600e-02j,  0.0000+0.0000e+00j,
                   0.0377-5.4600e-02j, -0.0698-1.0140e-01j,  0.5406+0.0000e+00j]],
        <BLANKLINE>
                [[ 0.4803+0.0000e+00j, -0.0923-1.2690e-01j,  0.1221+3.1900e-02j,
                   0.0000+0.0000e+00j, -0.0532-1.0000e-04j,  0.0830-8.8000e-03j],
                 [-0.0923+1.2690e-01j,  0.4942+0.0000e+00j,  0.0287+3.6500e-02j,
                   0.0532+1.0000e-04j,  0.0000+0.0000e+00j, -0.0995+4.8100e-02j],
                 [ 0.1221-3.1900e-02j,  0.0287-3.6500e-02j,  0.6290+0.0000e+00j,
                  -0.0830+8.8000e-03j,  0.0995-4.8100e-02j,  0.0000+0.0000e+00j],
                 [ 0.0000+0.0000e+00j,  0.0532-1.0000e-04j, -0.0830-8.8000e-03j,
                   0.4803+0.0000e+00j, -0.0923+1.2690e-01j,  0.1221-3.1900e-02j],
                 [-0.0532+1.0000e-04j,  0.0000+0.0000e+00j,  0.0995+4.8100e-02j,
                  -0.0923-1.2690e-01j,  0.4942+0.0000e+00j,  0.0287-3.6500e-02j],
                 [ 0.0830+8.8000e-03j, -0.0995-4.8100e-02j,  0.0000+0.0000e+00j,
                   0.1221+3.1900e-02j,  0.0287+3.6500e-02j,  0.6290+0.0000e+00j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if delta is None: 
        delta = torch.arange(1,n+1,device=device,dtype=torch.get_default_dtype())
        delta = delta/torch.linalg.norm(delta)
    assert delta.shape==(n,)
    assert not delta.is_complex()
    if rand_Spn is None: rand_Spn = rand_Spn_SVD
    assert rand_Spn in [rand_Spn_SVD]
    assert not torch.is_complex(delta)
    q = rand_Spn(N=N,n=n,k=n,seed=seed,device=delta.device,qp_unif_gen=qp_unif_gen)
    lam_paired = torch.cat([delta, delta], dim=-1)
    x = torch.einsum("...ij,...j,...kj->...ik", q,lam_paired.to(q.dtype),q.conj())
    return x

rand_LGr_R 

rand_LGr_R(N, n, seed=None, rand_Un=None, qp_unif_gen=None, device=None)

Generate a batch of N random real Lagrangian Grassmannian matrices of size (N, n, n).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Dimension of the Lagrangian subspace.

 required
seed int

Random seed for reproducibility.

 None
rand_Un callable

Function to generate random unitary matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random real Lagrangian Grassmannian matrices of size (N, n, n).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_LGr_R(2,3,seed=7)
>>> x.shape 
torch.Size([2, 3, 3])
>>> x
tensor([[[ 0.5347-0.6183j, -0.4216+0.1853j,  0.3367+0.0794j],
         [-0.4216+0.1853j, -0.7076-0.3250j,  0.3573-0.2321j],
         [ 0.3367+0.0794j,  0.3573-0.2321j,  0.3638-0.7526j]],

        [[ 0.1151+0.8999j, -0.0811+0.0124j,  0.3417+0.2311j],
         [-0.0811+0.0124j,  0.1517+0.9848j, -0.0134+0.0125j],
         [ 0.3417+0.2311j, -0.0134+0.0125j, -0.8817+0.2282j]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
def rand_LGr_R(N, n, seed=None, rand_Un=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random real Lagrangian Grassmannian matrices of size `(N, n, n)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Dimension of the Lagrangian subspace.
        seed (int): Random seed for reproducibility.
        rand_Un (callable): Function to generate random unitary matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random real Lagrangian Grassmannian matrices of size `(N, n, n)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_LGr_R(2,3,seed=7)
        >>> x.shape 
        torch.Size([2, 3, 3])
        >>> x
        tensor([[[ 0.5347-0.6183j, -0.4216+0.1853j,  0.3367+0.0794j],
                 [-0.4216+0.1853j, -0.7076-0.3250j,  0.3573-0.2321j],
                 [ 0.3367+0.0794j,  0.3573-0.2321j,  0.3638-0.7526j]],
        <BLANKLINE>
                [[ 0.1151+0.8999j, -0.0811+0.0124j,  0.3417+0.2311j],
                 [-0.0811+0.0124j,  0.1517+0.9848j, -0.0134+0.0125j],
                 [ 0.3417+0.2311j, -0.0134+0.0125j, -0.8817+0.2282j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if rand_Un is None: rand_Un = rand_Un_QR
    assert rand_Un in [rand_Un_QR,rand_Un_eig]
    q = rand_Un(N=N,n=n,k=n,seed=seed,device=device,qp_unif_gen=qp_unif_gen)
    x = torch.einsum("...ij,...kj->...ik",q,q)
    return x

rand_LGr_C 

rand_LGr_C(N, n, seed=None, rand_Spn=None, qp_unif_gen=None, device=None)

Generate a batch of N random complex Lagrangian Grassmannian matrices of size (N, 2n, 2n).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Half the dimension of the square matrix.

 required
seed int

Random seed for reproducibility.

 None
rand_Spn callable

Function to generate random symplectic unitary matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random complex Lagrangian Grassmannian matrices of size (N, 2n, 2n).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_LGr_C(2,3,seed=7)
>>> x.shape 
torch.Size([2, 6, 6])
>>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
tensor([[[ 0.2639+0.0000j,  0.0442+0.5590j,  0.4301+0.3442j, -0.3209-0.2230j,
          -0.3178+0.0323j, -0.2092-0.1182j],
         [ 0.0442-0.5590j, -0.0769+0.0000j, -0.2988+0.1905j, -0.3178+0.0323j,
          -0.2813-0.1818j,  0.3129+0.4919j],
         [ 0.4301-0.3442j, -0.2988-0.1905j, -0.0221+0.0000j, -0.2092-0.1182j,
           0.3129+0.4919j,  0.1832-0.3732j],
         [ 0.3209-0.2230j,  0.3178+0.0323j,  0.2092-0.1182j,  0.2639+0.0000j,
           0.0442-0.5590j,  0.4301-0.3442j],
         [ 0.3178+0.0323j,  0.2813-0.1818j, -0.3129+0.4919j,  0.0442+0.5590j,
          -0.0769+0.0000j, -0.2988-0.1905j],
         [ 0.2092-0.1182j, -0.3129+0.4919j, -0.1832-0.3732j,  0.4301+0.3442j,
          -0.2988+0.1905j, -0.0221+0.0000j]],

        [[-0.0856+0.0000j, -0.1772+0.0371j, -0.3202+0.0526j, -0.3455+0.5465j,
          -0.4689+0.1672j,  0.4316-0.0500j],
         [-0.1772-0.0371j, -0.0281+0.0000j,  0.0949-0.2171j, -0.4689+0.1672j,
           0.2348+0.5843j, -0.4891+0.1633j],
         [-0.3202-0.0526j,  0.0949+0.2171j,  0.0069+0.0000j,  0.4316-0.0500j,
          -0.4891+0.1633j, -0.1544+0.6000j],
         [ 0.3455+0.5465j,  0.4689+0.1672j, -0.4316-0.0500j, -0.0856+0.0000j,
          -0.1772-0.0371j, -0.3202-0.0526j],
         [ 0.4689+0.1672j, -0.2348+0.5843j,  0.4891+0.1633j, -0.1772+0.0371j,
          -0.0281+0.0000j,  0.0949+0.2171j],
         [-0.4316-0.0500j,  0.4891+0.1633j,  0.1544+0.6000j, -0.3202+0.0526j,
           0.0949-0.2171j,  0.0069+0.0000j]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
def rand_LGr_C(N, n, seed=None, rand_Spn=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random complex Lagrangian Grassmannian matrices of size `(N, 2n, 2n)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Half the dimension of the square matrix.
        seed (int): Random seed for reproducibility.
        rand_Spn (callable): Function to generate random symplectic unitary matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random complex Lagrangian Grassmannian matrices of size `(N, 2n, 2n)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_LGr_C(2,3,seed=7)
        >>> x.shape 
        torch.Size([2, 6, 6])
        >>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
        tensor([[[ 0.2639+0.0000j,  0.0442+0.5590j,  0.4301+0.3442j, -0.3209-0.2230j,
                  -0.3178+0.0323j, -0.2092-0.1182j],
                 [ 0.0442-0.5590j, -0.0769+0.0000j, -0.2988+0.1905j, -0.3178+0.0323j,
                  -0.2813-0.1818j,  0.3129+0.4919j],
                 [ 0.4301-0.3442j, -0.2988-0.1905j, -0.0221+0.0000j, -0.2092-0.1182j,
                   0.3129+0.4919j,  0.1832-0.3732j],
                 [ 0.3209-0.2230j,  0.3178+0.0323j,  0.2092-0.1182j,  0.2639+0.0000j,
                   0.0442-0.5590j,  0.4301-0.3442j],
                 [ 0.3178+0.0323j,  0.2813-0.1818j, -0.3129+0.4919j,  0.0442+0.5590j,
                  -0.0769+0.0000j, -0.2988-0.1905j],
                 [ 0.2092-0.1182j, -0.3129+0.4919j, -0.1832-0.3732j,  0.4301+0.3442j,
                  -0.2988+0.1905j, -0.0221+0.0000j]],
        <BLANKLINE>
                [[-0.0856+0.0000j, -0.1772+0.0371j, -0.3202+0.0526j, -0.3455+0.5465j,
                  -0.4689+0.1672j,  0.4316-0.0500j],
                 [-0.1772-0.0371j, -0.0281+0.0000j,  0.0949-0.2171j, -0.4689+0.1672j,
                   0.2348+0.5843j, -0.4891+0.1633j],
                 [-0.3202-0.0526j,  0.0949+0.2171j,  0.0069+0.0000j,  0.4316-0.0500j,
                  -0.4891+0.1633j, -0.1544+0.6000j],
                 [ 0.3455+0.5465j,  0.4689+0.1672j, -0.4316-0.0500j, -0.0856+0.0000j,
                  -0.1772-0.0371j, -0.3202-0.0526j],
                 [ 0.4689+0.1672j, -0.2348+0.5843j,  0.4891+0.1633j, -0.1772+0.0371j,
                  -0.0281+0.0000j,  0.0949+0.2171j],
                 [-0.4316-0.0500j,  0.4891+0.1633j,  0.1544+0.6000j, -0.3202+0.0526j,
                   0.0949-0.2171j,  0.0069+0.0000j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if rand_Spn is None: rand_Spn = rand_Spn_SVD
    assert rand_Spn in [rand_Spn_SVD]
    q = rand_Spn(N=N,n=n,k=n,seed=seed,device=device,qp_unif_gen=qp_unif_gen)
    I_nn = torch.diag(torch.cat([
        torch.ones(n,device=device), 
        -torch.ones(n,device=device)
    ])).to(q.dtype)
    q_left_corner = torch.einsum("ij,...kj,kl->...il",I_nn,q.conj(),I_nn)
    x = torch.einsum("...ij,...jk->...ik",q,q_left_corner)
    return x

rand_LGr_H 

rand_LGr_H(N, n, seed=None, rand_Un=None, qp_unif_gen=None, device=None)

Generate a batch of N random quaternionic Lagrangian Grassmannian matrices of size (N, 2n, 2n).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Half the dimension of the square matrix.

 required
seed int

Random seed for reproducibility.

 None
rand_Un callable

Function to generate random unitary matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random quaternionic Lagrangian Grassmannian matrices of size (N, 2n, 2n).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_LGr_H(2,3,seed=7)
>>> x.shape 
torch.Size([2, 6, 6])
>>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
tensor([[[ 0.2525+0.4835j,  0.3498+0.2373j,  0.0075+0.2605j,  0.0000+0.0000j,
          -0.0722+0.0447j, -0.0564+0.6675j],
         [-0.1068+0.5171j, -0.1195+0.0116j,  0.6605+0.3107j,  0.0722-0.0447j,
           0.0000+0.0000j,  0.1007-0.3958j],
         [-0.0877-0.3649j,  0.1217+0.4282j,  0.1840+0.1076j,  0.0564-0.6675j,
          -0.1007+0.3958j,  0.0000+0.0000j],
         [ 0.0000+0.0000j,  0.5042+0.1211j, -0.1175+0.0137j,  0.2525+0.4835j,
          -0.1068+0.5171j, -0.0877-0.3649j],
         [-0.5042-0.1211j,  0.0000+0.0000j,  0.4195-0.4047j,  0.3498+0.2373j,
          -0.1195+0.0116j,  0.1217+0.4282j],
         [ 0.1175-0.0137j, -0.4195+0.4047j,  0.0000+0.0000j,  0.0075+0.2605j,
           0.6605+0.3107j,  0.1840+0.1076j]],

        [[ 0.0364+0.3162j, -0.0919+0.2657j, -0.3184+0.4507j,  0.0000+0.0000j,
          -0.0340-0.2485j, -0.1330+0.6592j],
         [-0.2394+0.0866j,  0.0133-0.2684j, -0.6471+0.0804j,  0.0340+0.2485j,
           0.0000+0.0000j,  0.5890-0.1670j],
         [ 0.0843-0.2270j, -0.1381-0.0200j, -0.3063-0.0327j,  0.1330-0.6592j,
          -0.5890+0.1670j,  0.0000+0.0000j],
         [ 0.0000+0.0000j, -0.1605+0.8300j, -0.1159-0.2173j,  0.0364+0.3162j,
          -0.2394+0.0866j,  0.0843-0.2270j],
         [ 0.1605-0.8300j,  0.0000+0.0000j, -0.1564+0.3004j, -0.0919+0.2657j,
           0.0133-0.2684j, -0.1381-0.0200j],
         [ 0.1159+0.2173j,  0.1564-0.3004j,  0.0000+0.0000j, -0.3184+0.4507j,
          -0.6471+0.0804j, -0.3063-0.0327j]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
def rand_LGr_H(N, n, seed=None, rand_Un=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random quaternionic Lagrangian Grassmannian matrices of size `(N, 2n, 2n)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Half the dimension of the square matrix.
        seed (int): Random seed for reproducibility.
        rand_Un (callable): Function to generate random unitary matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random quaternionic Lagrangian Grassmannian matrices of size `(N, 2n, 2n)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_LGr_H(2,3,seed=7)
        >>> x.shape 
        torch.Size([2, 6, 6])
        >>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
        tensor([[[ 0.2525+0.4835j,  0.3498+0.2373j,  0.0075+0.2605j,  0.0000+0.0000j,
                  -0.0722+0.0447j, -0.0564+0.6675j],
                 [-0.1068+0.5171j, -0.1195+0.0116j,  0.6605+0.3107j,  0.0722-0.0447j,
                   0.0000+0.0000j,  0.1007-0.3958j],
                 [-0.0877-0.3649j,  0.1217+0.4282j,  0.1840+0.1076j,  0.0564-0.6675j,
                  -0.1007+0.3958j,  0.0000+0.0000j],
                 [ 0.0000+0.0000j,  0.5042+0.1211j, -0.1175+0.0137j,  0.2525+0.4835j,
                  -0.1068+0.5171j, -0.0877-0.3649j],
                 [-0.5042-0.1211j,  0.0000+0.0000j,  0.4195-0.4047j,  0.3498+0.2373j,
                  -0.1195+0.0116j,  0.1217+0.4282j],
                 [ 0.1175-0.0137j, -0.4195+0.4047j,  0.0000+0.0000j,  0.0075+0.2605j,
                   0.6605+0.3107j,  0.1840+0.1076j]],
        <BLANKLINE>
                [[ 0.0364+0.3162j, -0.0919+0.2657j, -0.3184+0.4507j,  0.0000+0.0000j,
                  -0.0340-0.2485j, -0.1330+0.6592j],
                 [-0.2394+0.0866j,  0.0133-0.2684j, -0.6471+0.0804j,  0.0340+0.2485j,
                   0.0000+0.0000j,  0.5890-0.1670j],
                 [ 0.0843-0.2270j, -0.1381-0.0200j, -0.3063-0.0327j,  0.1330-0.6592j,
                  -0.5890+0.1670j,  0.0000+0.0000j],
                 [ 0.0000+0.0000j, -0.1605+0.8300j, -0.1159-0.2173j,  0.0364+0.3162j,
                  -0.2394+0.0866j,  0.0843-0.2270j],
                 [ 0.1605-0.8300j,  0.0000+0.0000j, -0.1564+0.3004j, -0.0919+0.2657j,
                   0.0133-0.2684j, -0.1381-0.0200j],
                 [ 0.1159+0.2173j,  0.1564-0.3004j,  0.0000+0.0000j, -0.3184+0.4507j,
                  -0.6471+0.0804j, -0.3063-0.0327j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if rand_Un is None: rand_Un = rand_Un_QR
    assert rand_Un in [rand_Un_QR,rand_Un_eig]
    q = rand_Un(N=N,n=2*n,k=2*n,seed=seed,device=device,qp_unif_gen=qp_unif_gen)
    I_N = torch.eye(n,device=device,dtype=q.dtype)
    O_N = torch.zeros(n,n,device=device,dtype=q.dtype)
    J = torch.cat([
        torch.cat([O_N,-I_N],dim=-1),
        torch.cat([I_N, O_N],dim=-1)],dim=-2)
    q_S = -torch.einsum("ij,...kj,kl->...il",J,q,J)
    x = torch.einsum("...ij,...jk->...ik",q,q_S)
    return x

rand_Stiefel_R 

rand_Stiefel_R(N, n, k=1, seed=None, rand_On=None, qp_unif_gen=None, device=None)

Generate a batch of N random real Stiefel matrices of size (N, n, k).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Number of rows.

 required
k int

Number of columns.

 1
seed int

Random seed for reproducibility.

 None
rand_On callable

Function to generate random orthogonal matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random real Stiefel matrices of size (N, n, k).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_Stiefel_R(2,4,3,seed=7)
>>> x.shape 
torch.Size([2, 4, 3])
>>> x
tensor([[[-0.7177,  0.3121,  0.1379],
         [ 0.1644,  0.3183, -0.9033],
         [ 0.1398, -0.7875, -0.0974],
         [-0.6620, -0.4256, -0.3944]],

        [[-0.7589,  0.1929,  0.3816],
         [ 0.1093, -0.6896,  0.7074],
         [ 0.6078,  0.5222,  0.4722],
         [ 0.2064, -0.4631, -0.3620]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
def rand_Stiefel_R(N, n, k=1, seed=None, rand_On=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random real Stiefel matrices of size `(N, n, k)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Number of rows.
        k (int): Number of columns.
        seed (int): Random seed for reproducibility.
        rand_On (callable): Function to generate random orthogonal matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random real Stiefel matrices of size `(N, n, k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_Stiefel_R(2,4,3,seed=7)
        >>> x.shape 
        torch.Size([2, 4, 3])
        >>> x
        tensor([[[-0.7177,  0.3121,  0.1379],
                 [ 0.1644,  0.3183, -0.9033],
                 [ 0.1398, -0.7875, -0.0974],
                 [-0.6620, -0.4256, -0.3944]],
        <BLANKLINE>
                [[-0.7589,  0.1929,  0.3816],
                 [ 0.1093, -0.6896,  0.7074],
                 [ 0.6078,  0.5222,  0.4722],
                 [ 0.2064, -0.4631, -0.3620]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if rand_On is None: rand_On = rand_On_QR
    assert rand_On in [rand_On_QR,rand_On_eig]
    assert 1<=k<=n
    x = rand_On(N=N,n=n,k=k,seed=seed,device=device,qp_unif_gen=qp_unif_gen)
    return x

rand_Stiefel_C 

rand_Stiefel_C(N, n, k=1, seed=None, rand_Un=None, qp_unif_gen=None, device=None)

Generate a batch of N random complex Stiefel matrices of size (N, n, k).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Number of rows.

 required
k int

Number of columns.

 1
seed int

Random seed for reproducibility.

 None
rand_Un callable

Function to generate random unitary matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random complex Stiefel matrices of size (N, n, k).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_Stiefel_C(2,4,3,seed=7)
>>> x.shape 
torch.Size([2, 4, 3])
>>> x
tensor([[[-0.6106+0.0789j,  0.0853-0.1095j,  0.2571-0.3622j],
         [ 0.1189-0.3923j,  0.2634-0.6302j, -0.1178-0.5404j],
         [-0.3917+0.2437j,  0.6549-0.2208j, -0.2803+0.4433j],
         [ 0.3137+0.3765j,  0.1450+0.1248j, -0.3241-0.3415j]],

        [[ 0.0893+0.1079j, -0.1481-0.1997j,  0.8927-0.1751j],
         [ 0.0739-0.3969j, -0.1797-0.3939j,  0.0373-0.1982j],
         [-0.0191-0.4580j,  0.0510-0.6749j, -0.1948+0.0873j],
         [ 0.3224-0.7095j, -0.1439+0.5215j,  0.1805+0.2314j]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
def rand_Stiefel_C(N, n, k=1, seed=None, rand_Un=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random complex Stiefel matrices of size `(N, n, k)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Number of rows.
        k (int): Number of columns.
        seed (int): Random seed for reproducibility.
        rand_Un (callable): Function to generate random unitary matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random complex Stiefel matrices of size `(N, n, k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_Stiefel_C(2,4,3,seed=7)
        >>> x.shape 
        torch.Size([2, 4, 3])
        >>> x
        tensor([[[-0.6106+0.0789j,  0.0853-0.1095j,  0.2571-0.3622j],
                 [ 0.1189-0.3923j,  0.2634-0.6302j, -0.1178-0.5404j],
                 [-0.3917+0.2437j,  0.6549-0.2208j, -0.2803+0.4433j],
                 [ 0.3137+0.3765j,  0.1450+0.1248j, -0.3241-0.3415j]],
        <BLANKLINE>
                [[ 0.0893+0.1079j, -0.1481-0.1997j,  0.8927-0.1751j],
                 [ 0.0739-0.3969j, -0.1797-0.3939j,  0.0373-0.1982j],
                 [-0.0191-0.4580j,  0.0510-0.6749j, -0.1948+0.0873j],
                 [ 0.3224-0.7095j, -0.1439+0.5215j,  0.1805+0.2314j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if rand_Un is None: rand_Un = rand_Un_QR
    assert rand_Un in [rand_Un_QR,rand_Un_eig]
    assert 1<=k<=n
    x = rand_Un(N=N,n=n,k=k,seed=seed,device=device,qp_unif_gen=qp_unif_gen)
    return x

rand_Stiefel_H 

rand_Stiefel_H(N, n, k=1, seed=None, rand_Spn=None, qp_unif_gen=None, device=None)

Generate a batch of N random quaternionic Stiefel matrices of size (N, 2n, 2k).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Half the number of rows.

 required
k int

Half the number of columns.

 1
seed int

Random seed for reproducibility.

 None
rand_Spn callable

Function to generate random symplectic unitary matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random quaternionic Stiefel matrices of size (N, 2n, 2k).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_Stiefel_H(2,4,3,seed=7)
>>> x.shape 
torch.Size([2, 8, 6])
>>> x
tensor([[[-0.4140+0.3030j,  0.0601+0.0313j,  0.2725-0.3912j,  0.0845+0.0474j,
           0.0467-0.2241j, -0.2175-0.1204j],
         [-0.5555+0.0748j,  0.3572+0.0449j, -0.0927+0.1001j,  0.1535+0.3413j,
          -0.2548-0.3059j,  0.1301+0.1163j],
         [-0.0531-0.1221j, -0.0396-0.3151j, -0.1012+0.2916j,  0.1274-0.2141j,
          -0.0340-0.2394j,  0.5060-0.4238j],
         [-0.3121-0.2754j, -0.4287-0.3945j, -0.0973-0.0722j, -0.0424-0.1354j,
          -0.3801-0.1076j, -0.1164+0.3205j],
         [-0.0845+0.0474j, -0.0467-0.2241j,  0.2175-0.1204j, -0.4140-0.3030j,
           0.0601-0.0313j,  0.2725+0.3912j],
         [-0.1535+0.3413j,  0.2548-0.3059j, -0.1301+0.1163j, -0.5555-0.0748j,
           0.3572-0.0449j, -0.0927-0.1001j],
         [-0.1274-0.2141j,  0.0340-0.2394j, -0.5060-0.4238j, -0.0531+0.1221j,
          -0.0396+0.3151j, -0.1012-0.2916j],
         [ 0.0424-0.1354j,  0.3801-0.1076j,  0.1164+0.3205j, -0.3121+0.2754j,
          -0.4287+0.3945j, -0.0973+0.0722j]],

        [[ 0.2331-0.3342j,  0.4754+0.1704j, -0.0965+0.1623j, -0.3163-0.3840j,
           0.2357+0.3150j, -0.2709+0.1039j],
         [ 0.1156-0.2160j,  0.0750-0.2084j, -0.0375+0.0934j, -0.0851-0.1242j,
          -0.1064-0.1801j,  0.3919+0.0141j],
         [ 0.3396+0.1853j,  0.3284-0.2578j, -0.0459+0.3762j,  0.2561+0.2939j,
          -0.3221-0.2715j, -0.3632+0.0764j],
         [ 0.2941-0.0642j,  0.1781+0.0226j,  0.4171+0.1087j, -0.3197-0.0961j,
          -0.2625-0.2109j,  0.1682-0.4698j],
         [ 0.3163-0.3840j, -0.2357+0.3150j,  0.2709+0.1039j,  0.2331+0.3342j,
           0.4754-0.1704j, -0.0965-0.1623j],
         [ 0.0851-0.1242j,  0.1064-0.1801j, -0.3919+0.0141j,  0.1156+0.2160j,
           0.0750+0.2084j, -0.0375-0.0934j],
         [-0.2561+0.2939j,  0.3221-0.2715j,  0.3632+0.0764j,  0.3396-0.1853j,
           0.3284+0.2578j, -0.0459-0.3762j],
         [ 0.3197-0.0961j,  0.2625-0.2109j, -0.1682-0.4698j,  0.2941+0.0642j,
           0.1781-0.0226j,  0.4171-0.1087j]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
def rand_Stiefel_H(N, n, k=1, seed=None, rand_Spn=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random quaternionic Stiefel matrices of size `(N, 2n, 2k)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Half the number of rows.
        k (int): Half the number of columns.
        seed (int): Random seed for reproducibility.
        rand_Spn (callable): Function to generate random symplectic unitary matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random quaternionic Stiefel matrices of size `(N, 2n, 2k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_Stiefel_H(2,4,3,seed=7)
        >>> x.shape 
        torch.Size([2, 8, 6])
        >>> x
        tensor([[[-0.4140+0.3030j,  0.0601+0.0313j,  0.2725-0.3912j,  0.0845+0.0474j,
                   0.0467-0.2241j, -0.2175-0.1204j],
                 [-0.5555+0.0748j,  0.3572+0.0449j, -0.0927+0.1001j,  0.1535+0.3413j,
                  -0.2548-0.3059j,  0.1301+0.1163j],
                 [-0.0531-0.1221j, -0.0396-0.3151j, -0.1012+0.2916j,  0.1274-0.2141j,
                  -0.0340-0.2394j,  0.5060-0.4238j],
                 [-0.3121-0.2754j, -0.4287-0.3945j, -0.0973-0.0722j, -0.0424-0.1354j,
                  -0.3801-0.1076j, -0.1164+0.3205j],
                 [-0.0845+0.0474j, -0.0467-0.2241j,  0.2175-0.1204j, -0.4140-0.3030j,
                   0.0601-0.0313j,  0.2725+0.3912j],
                 [-0.1535+0.3413j,  0.2548-0.3059j, -0.1301+0.1163j, -0.5555-0.0748j,
                   0.3572-0.0449j, -0.0927-0.1001j],
                 [-0.1274-0.2141j,  0.0340-0.2394j, -0.5060-0.4238j, -0.0531+0.1221j,
                  -0.0396+0.3151j, -0.1012-0.2916j],
                 [ 0.0424-0.1354j,  0.3801-0.1076j,  0.1164+0.3205j, -0.3121+0.2754j,
                  -0.4287+0.3945j, -0.0973+0.0722j]],
        <BLANKLINE>
                [[ 0.2331-0.3342j,  0.4754+0.1704j, -0.0965+0.1623j, -0.3163-0.3840j,
                   0.2357+0.3150j, -0.2709+0.1039j],
                 [ 0.1156-0.2160j,  0.0750-0.2084j, -0.0375+0.0934j, -0.0851-0.1242j,
                  -0.1064-0.1801j,  0.3919+0.0141j],
                 [ 0.3396+0.1853j,  0.3284-0.2578j, -0.0459+0.3762j,  0.2561+0.2939j,
                  -0.3221-0.2715j, -0.3632+0.0764j],
                 [ 0.2941-0.0642j,  0.1781+0.0226j,  0.4171+0.1087j, -0.3197-0.0961j,
                  -0.2625-0.2109j,  0.1682-0.4698j],
                 [ 0.3163-0.3840j, -0.2357+0.3150j,  0.2709+0.1039j,  0.2331+0.3342j,
                   0.4754-0.1704j, -0.0965-0.1623j],
                 [ 0.0851-0.1242j,  0.1064-0.1801j, -0.3919+0.0141j,  0.1156+0.2160j,
                   0.0750+0.2084j, -0.0375-0.0934j],
                 [-0.2561+0.2939j,  0.3221-0.2715j,  0.3632+0.0764j,  0.3396-0.1853j,
                   0.3284+0.2578j, -0.0459-0.3762j],
                 [ 0.3197-0.0961j,  0.2625-0.2109j, -0.1682-0.4698j,  0.2941+0.0642j,
                   0.1781-0.0226j,  0.4171-0.1087j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if rand_Spn is None: rand_Spn = rand_Spn_SVD
    assert rand_Spn in [rand_Spn_SVD]
    assert 1<=k<=n
    x = rand_Spn(N=N,n=n,k=k,seed=seed,device=device,qp_unif_gen=qp_unif_gen)
    return x

rand_Gr_R 

rand_Gr_R(N, n, k=1, seed=None, rand_On=None, qp_unif_gen=None, device=None)

Generate a batch of N random real Grassmannian matrices of size (N, n, n).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Dimension of the space.

 required
k int

Dimension of the subspace.

 1
seed int

Random seed for reproducibility.

 None
rand_On callable

Function to generate random orthogonal matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random real Grassmannian matrices of size (N, n, n).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_Gr_R(2,4,3,seed=7)
>>> x.shape 
torch.Size([2, 4, 4])
>>> x
tensor([[[ 0.2631, -0.2864, -0.7190,  0.5760],
         [-0.2864,  0.8887, -0.2795,  0.2239],
         [-0.7190, -0.2795,  0.2984,  0.5620],
         [ 0.5760,  0.2239,  0.5620,  0.5498]],

        [[ 0.5176,  0.1079, -0.3608, -0.7683],
         [ 0.1079,  0.9759,  0.0807,  0.1718],
         [-0.3608,  0.0807,  0.7302, -0.5746],
         [-0.7683,  0.1718, -0.5746, -0.2237]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
def rand_Gr_R(N, n, k=1, seed=None, rand_On=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random real Grassmannian matrices of size `(N, n, n)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Dimension of the space.
        k (int): Dimension of the subspace.
        seed (int): Random seed for reproducibility.
        rand_On (callable): Function to generate random orthogonal matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random real Grassmannian matrices of size `(N, n, n)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_Gr_R(2,4,3,seed=7)
        >>> x.shape 
        torch.Size([2, 4, 4])
        >>> x
        tensor([[[ 0.2631, -0.2864, -0.7190,  0.5760],
                 [-0.2864,  0.8887, -0.2795,  0.2239],
                 [-0.7190, -0.2795,  0.2984,  0.5620],
                 [ 0.5760,  0.2239,  0.5620,  0.5498]],
        <BLANKLINE>
                [[ 0.5176,  0.1079, -0.3608, -0.7683],
                 [ 0.1079,  0.9759,  0.0807,  0.1718],
                 [-0.3608,  0.0807,  0.7302, -0.5746],
                 [-0.7683,  0.1718, -0.5746, -0.2237]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if rand_On is None: rand_On = rand_On_QR
    assert rand_On in [rand_On_QR,rand_On_eig]
    assert 1<=k<=n
    q = rand_On(N=N,n=n,k=k,seed=seed,qp_unif_gen=qp_unif_gen)
    x = 2*torch.einsum("...ij,...kj->...ik",q,q)-torch.eye(n)
    return x

rand_Gr_C 

rand_Gr_C(N, n, k=1, seed=None, rand_Un=None, qp_unif_gen=None, device=None)

Generate a batch of N random complex Grassmannian matrices of size (N, n, n).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Dimension of the space.

 required
k int

Dimension of the subspace.

 1
seed int

Random seed for reproducibility.

 None
rand_Un callable

Function to generate random unitary matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random complex Grassmannian matrices of size (N, n, n).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_Gr_C(2,4,3,seed=7)
>>> x.shape 
torch.Size([2, 4, 4])
>>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
tensor([[[ 0.1911+0.0000j,  0.3068-0.0473j,  0.2115+0.1051j, -0.2456+0.8666j],
         [ 0.3068+0.0473j,  0.8809+0.0000j, -0.0741-0.0522j,  0.1438-0.3143j],
         [ 0.2115-0.1051j, -0.0741+0.0522j,  0.9310+0.0000j, -0.0484-0.2585j],
         [-0.2456-0.8666j,  0.1438+0.3143j, -0.0484+0.2585j, -0.0031+0.0000j]],

        [[ 0.8180+0.0000j,  0.2741+0.3828j, -0.2263-0.2302j, -0.0199-0.0682j],
         [ 0.2741-0.3828j, -0.2178+0.0000j,  0.8249-0.1292j,  0.1734+0.0610j],
         [-0.2263+0.2302j,  0.8249+0.1292j,  0.4275+0.0000j, -0.1110-0.0597j],
         [-0.0199+0.0682j,  0.1734-0.0610j, -0.1110+0.0597j,  0.9723+0.0000j]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
def rand_Gr_C(N, n, k=1, seed=None, rand_Un=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random complex Grassmannian matrices of size `(N, n, n)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Dimension of the space.
        k (int): Dimension of the subspace.
        seed (int): Random seed for reproducibility.
        rand_Un (callable): Function to generate random unitary matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random complex Grassmannian matrices of size `(N, n, n)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_Gr_C(2,4,3,seed=7)
        >>> x.shape 
        torch.Size([2, 4, 4])
        >>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
        tensor([[[ 0.1911+0.0000j,  0.3068-0.0473j,  0.2115+0.1051j, -0.2456+0.8666j],
                 [ 0.3068+0.0473j,  0.8809+0.0000j, -0.0741-0.0522j,  0.1438-0.3143j],
                 [ 0.2115-0.1051j, -0.0741+0.0522j,  0.9310+0.0000j, -0.0484-0.2585j],
                 [-0.2456-0.8666j,  0.1438+0.3143j, -0.0484+0.2585j, -0.0031+0.0000j]],
        <BLANKLINE>
                [[ 0.8180+0.0000j,  0.2741+0.3828j, -0.2263-0.2302j, -0.0199-0.0682j],
                 [ 0.2741-0.3828j, -0.2178+0.0000j,  0.8249-0.1292j,  0.1734+0.0610j],
                 [-0.2263+0.2302j,  0.8249+0.1292j,  0.4275+0.0000j, -0.1110-0.0597j],
                 [-0.0199+0.0682j,  0.1734-0.0610j, -0.1110+0.0597j,  0.9723+0.0000j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if rand_Un is None: rand_Un = rand_Un_QR
    assert rand_Un in [rand_Un_QR,rand_Un_eig]
    assert 1<=k<=n
    q = rand_Un(N=N,n=n,k=k,seed=seed,qp_unif_gen=qp_unif_gen)
    x = 2*torch.einsum("...ij,...kj->...ik",q,q.conj())-torch.eye(n)
    return x

rand_Gr_H 

rand_Gr_H(N, n, k=1, seed=None, rand_Spn=None, qp_unif_gen=None, device=None)

Generate a batch of N random quaternionic Grassmannian matrices of size (N, 2n, 2n).

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Half the dimension of the space.

 required
k int

Half the dimension of the subspace.

 1
seed int

Random seed for reproducibility.

 None
rand_Spn callable

Function to generate random symplectic unitary matrices.

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None
device str

Device to store the tensor on.

 None

Returns:

Name Type Description
x Tensor

A batch of N random quaternionic Grassmannian matrices of size (N, 2n, 2n).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> x = rand_Gr_H(2,4,3,seed=7)
>>> x.shape 
torch.Size([2, 8, 8])
>>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
tensor([[[ 0.2373+0.0000j,  0.5093-0.1209j, -0.3505-0.3980j, -0.0150+0.0857j,
           0.0000+0.0000j,  0.2004+0.0635j, -0.0076+0.2635j, -0.4841-0.1813j],
         [ 0.5093+0.1209j,  0.5828+0.0000j,  0.1510+0.3912j,  0.1659-0.0622j,
          -0.2004-0.0635j,  0.0000+0.0000j, -0.0784-0.0994j,  0.2977+0.1740j],
         [-0.3505+0.3980j,  0.1510-0.3912j,  0.5401+0.0000j,  0.0957+0.2163j,
           0.0076-0.2635j,  0.0784+0.0994j,  0.0000+0.0000j, -0.2873+0.1649j],
         [-0.0150-0.0857j,  0.1659+0.0622j,  0.0957-0.2163j,  0.6397+0.0000j,
           0.4841+0.1813j, -0.2977-0.1740j,  0.2873-0.1649j,  0.0000+0.0000j],
         [ 0.0000+0.0000j, -0.2004+0.0635j,  0.0076+0.2635j,  0.4841-0.1813j,
           0.2373+0.0000j,  0.5093+0.1209j, -0.3505+0.3980j, -0.0150-0.0857j],
         [ 0.2004-0.0635j,  0.0000+0.0000j,  0.0784-0.0994j, -0.2977+0.1740j,
           0.5093-0.1209j,  0.5828+0.0000j,  0.1510-0.3912j,  0.1659+0.0622j],
         [-0.0076-0.2635j, -0.0784+0.0994j,  0.0000+0.0000j,  0.2873+0.1649j,
          -0.3505-0.3980j,  0.1510+0.3912j,  0.5401+0.0000j,  0.0957-0.2163j],
         [-0.4841+0.1813j,  0.2977-0.1740j, -0.2873-0.1649j,  0.0000+0.0000j,
          -0.0150+0.0857j,  0.1659-0.0622j,  0.0957+0.2163j,  0.6397+0.0000j]],

        [[ 0.8864+0.0000j,  0.0124+0.3467j, -0.1082-0.0184j,  0.1425-0.0718j,
           0.0000+0.0000j,  0.1712+0.0225j,  0.0439+0.0287j, -0.1433-0.0589j],
         [ 0.0124-0.3467j, -0.3211+0.0000j, -0.0039-0.2936j,  0.4309+0.3820j,
          -0.1712-0.0225j,  0.0000+0.0000j, -0.2588+0.1369j,  0.3956-0.2944j],
         [-0.1082+0.0184j, -0.0039+0.2936j,  0.8699+0.0000j,  0.1942-0.0780j,
          -0.0439-0.0287j,  0.2588-0.1369j,  0.0000+0.0000j, -0.1091+0.0307j],
         [ 0.1425+0.0718j,  0.4309-0.3820j,  0.1942+0.0780j,  0.5649+0.0000j,
           0.1433+0.0589j, -0.3956+0.2944j,  0.1091-0.0307j,  0.0000+0.0000j],
         [ 0.0000+0.0000j, -0.1712+0.0225j, -0.0439+0.0287j,  0.1433-0.0589j,
           0.8864+0.0000j,  0.0124-0.3467j, -0.1082+0.0184j,  0.1425+0.0718j],
         [ 0.1712-0.0225j,  0.0000+0.0000j,  0.2588+0.1369j, -0.3956-0.2944j,
           0.0124+0.3467j, -0.3211+0.0000j, -0.0039+0.2936j,  0.4309-0.3820j],
         [ 0.0439-0.0287j, -0.2588-0.1369j,  0.0000+0.0000j,  0.1091+0.0307j,
          -0.1082-0.0184j, -0.0039-0.2936j,  0.8699+0.0000j,  0.1942+0.0780j],
         [-0.1433+0.0589j,  0.3956+0.2944j, -0.1091-0.0307j,  0.0000+0.0000j,
           0.1425-0.0718j,  0.4309+0.3820j,  0.1942-0.0780j,  0.5649+0.0000j]]])
Source code in rmatqmcpy/rand_Stiefel_Gr_flag_LGr.py 
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
def rand_Gr_H(N, n, k=1, seed=None, rand_Spn=None, qp_unif_gen=None, device=None):
    r"""
    Generate a batch of `N` random quaternionic Grassmannian matrices of size `(N, 2n, 2n)`.

    Args:
        N (int): Number of samples to generate.
        n (int): Half the dimension of the space.
        k (int): Half the dimension of the subspace.
        seed (int): Random seed for reproducibility.
        rand_Spn (callable): Function to generate random symplectic unitary matrices.
        qp_unif_gen (qmcpy.DiscreteDistribution): QMCPy distribution generator.
        device (str): Device to store the tensor on.

    Returns:
        x (torch.Tensor): A batch of `N` random quaternionic Grassmannian matrices of size `(N, 2n, 2n)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> x = rand_Gr_H(2,4,3,seed=7)
        >>> x.shape 
        torch.Size([2, 8, 8])
        >>> torch.complex(x.real.round(decimals=4)+0.,x.imag.round(decimals=4)+0.)
        tensor([[[ 0.2373+0.0000j,  0.5093-0.1209j, -0.3505-0.3980j, -0.0150+0.0857j,
                   0.0000+0.0000j,  0.2004+0.0635j, -0.0076+0.2635j, -0.4841-0.1813j],
                 [ 0.5093+0.1209j,  0.5828+0.0000j,  0.1510+0.3912j,  0.1659-0.0622j,
                  -0.2004-0.0635j,  0.0000+0.0000j, -0.0784-0.0994j,  0.2977+0.1740j],
                 [-0.3505+0.3980j,  0.1510-0.3912j,  0.5401+0.0000j,  0.0957+0.2163j,
                   0.0076-0.2635j,  0.0784+0.0994j,  0.0000+0.0000j, -0.2873+0.1649j],
                 [-0.0150-0.0857j,  0.1659+0.0622j,  0.0957-0.2163j,  0.6397+0.0000j,
                   0.4841+0.1813j, -0.2977-0.1740j,  0.2873-0.1649j,  0.0000+0.0000j],
                 [ 0.0000+0.0000j, -0.2004+0.0635j,  0.0076+0.2635j,  0.4841-0.1813j,
                   0.2373+0.0000j,  0.5093+0.1209j, -0.3505+0.3980j, -0.0150-0.0857j],
                 [ 0.2004-0.0635j,  0.0000+0.0000j,  0.0784-0.0994j, -0.2977+0.1740j,
                   0.5093-0.1209j,  0.5828+0.0000j,  0.1510-0.3912j,  0.1659+0.0622j],
                 [-0.0076-0.2635j, -0.0784+0.0994j,  0.0000+0.0000j,  0.2873+0.1649j,
                  -0.3505-0.3980j,  0.1510+0.3912j,  0.5401+0.0000j,  0.0957-0.2163j],
                 [-0.4841+0.1813j,  0.2977-0.1740j, -0.2873-0.1649j,  0.0000+0.0000j,
                  -0.0150+0.0857j,  0.1659-0.0622j,  0.0957+0.2163j,  0.6397+0.0000j]],
        <BLANKLINE>
                [[ 0.8864+0.0000j,  0.0124+0.3467j, -0.1082-0.0184j,  0.1425-0.0718j,
                   0.0000+0.0000j,  0.1712+0.0225j,  0.0439+0.0287j, -0.1433-0.0589j],
                 [ 0.0124-0.3467j, -0.3211+0.0000j, -0.0039-0.2936j,  0.4309+0.3820j,
                  -0.1712-0.0225j,  0.0000+0.0000j, -0.2588+0.1369j,  0.3956-0.2944j],
                 [-0.1082+0.0184j, -0.0039+0.2936j,  0.8699+0.0000j,  0.1942-0.0780j,
                  -0.0439-0.0287j,  0.2588-0.1369j,  0.0000+0.0000j, -0.1091+0.0307j],
                 [ 0.1425+0.0718j,  0.4309-0.3820j,  0.1942+0.0780j,  0.5649+0.0000j,
                   0.1433+0.0589j, -0.3956+0.2944j,  0.1091-0.0307j,  0.0000+0.0000j],
                 [ 0.0000+0.0000j, -0.1712+0.0225j, -0.0439+0.0287j,  0.1433-0.0589j,
                   0.8864+0.0000j,  0.0124-0.3467j, -0.1082+0.0184j,  0.1425+0.0718j],
                 [ 0.1712-0.0225j,  0.0000+0.0000j,  0.2588+0.1369j, -0.3956-0.2944j,
                   0.0124+0.3467j, -0.3211+0.0000j, -0.0039+0.2936j,  0.4309-0.3820j],
                 [ 0.0439-0.0287j, -0.2588-0.1369j,  0.0000+0.0000j,  0.1091+0.0307j,
                  -0.1082-0.0184j, -0.0039-0.2936j,  0.8699+0.0000j,  0.1942+0.0780j],
                 [-0.1433+0.0589j,  0.3956+0.2944j, -0.1091-0.0307j,  0.0000+0.0000j,
                   0.1425-0.0718j,  0.4309+0.3820j,  0.1942-0.0780j,  0.5649+0.0000j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    if rand_Spn is None: rand_Spn = rand_Spn_SVD
    assert rand_Spn in [rand_Spn_SVD]
    assert 1<=k<=n
    q = rand_Spn(N=N,n=n,k=k,seed=seed,qp_unif_gen=qp_unif_gen)
    x = 2*torch.einsum("...ij,...kj->...ik",q,q.conj())-torch.eye(2*n)
    return x

Random \(\mathrm{O}(n)\), \(\mathrm{U}(n)\) and \(\mathrm{Sp}(n)\) matrices on \(\mathbb{R}\), \(\mathbb{C}\), and \(\mathbb{H}\) respectively

rand_On_QR 

rand_On_QR(N, n, k, seed=None, device=None, qp_unif_gen=None)

Generate a batch of N random orthogonal matrices of size (N, n, k) using the QR decomposition.

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Number of rows.

 required
k int

Number of columns with k<=n.

 required
seed int

Random seed for reproducibility.

 None
device str

Device to store the tensor on (e.g., "cpu", "cuda").

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None

Returns:

Name Type Description
x Tensor

A batch of N random orthogonal matrices of size (N, n, k).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> q = rand_On_QR(2,3,3,seed=7)
>>> q.shape 
torch.Size([2, 3, 3])
>>> q
tensor([[[-0.9577, -0.0582, -0.2820],
         [ 0.2193,  0.4869, -0.8455],
         [ 0.1865, -0.8715, -0.4535]],

        [[-0.8182, -0.3566, -0.4510],
         [-0.5690,  0.3904,  0.7237],
         [ 0.0820, -0.8488,  0.5224]]])

>>> q = rand_On_QR(2,3,3,seed=7,qp_unif_gen=qp.Net)
>>> q.shape 
torch.Size([2, 3, 3])
>>> q
tensor([[[-0.1323, -0.4299,  0.8931],
         [-0.6249,  0.7356,  0.2615],
         [-0.7694, -0.5236, -0.3660]],

        [[ 0.8600, -0.4448, -0.2502],
         [ 0.4598,  0.8880,  0.0014],
         [ 0.2215, -0.1162,  0.9682]]])
Source code in rmatqmcpy/rand_On_Un_Spn.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def rand_On_QR(N, n, k, seed=None, device=None, qp_unif_gen=None):
    r"""
    Generate a batch of `N` random orthogonal matrices of size `(N, n, k)` using the QR decomposition. 

    Args:
        N (int): Number of samples to generate.
        n (int): Number of rows.
        k (int): Number of columns with `k<=n`.
        seed (int, optional): Random seed for reproducibility.
        device (str, optional): Device to store the tensor on (e.g., "cpu", "cuda").
        qp_unif_gen (qmcpy.DiscreteDistribution, optional): QMCPy distribution generator.

    Returns:
        x (torch.Tensor): A batch of `N` random orthogonal matrices of size `(N, n, k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 


        >>> q = rand_On_QR(2,3,3,seed=7)
        >>> q.shape 
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[-0.9577, -0.0582, -0.2820],
                 [ 0.2193,  0.4869, -0.8455],
                 [ 0.1865, -0.8715, -0.4535]],
        <BLANKLINE>
                [[-0.8182, -0.3566, -0.4510],
                 [-0.5690,  0.3904,  0.7237],
                 [ 0.0820, -0.8488,  0.5224]]])

        >>> q = rand_On_QR(2,3,3,seed=7,qp_unif_gen=qp.Net)
        >>> q.shape 
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[-0.1323, -0.4299,  0.8931],
                 [-0.6249,  0.7356,  0.2615],
                 [-0.7694, -0.5236, -0.3660]],
        <BLANKLINE>
                [[ 0.8600, -0.4448, -0.2502],
                 [ 0.4598,  0.8880,  0.0014],
                 [ 0.2215, -0.1162,  0.9682]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    rng = agsutil.get_torch_rng(seed,device=device)
    u = torch.from_numpy(qp_unif_gen(dimension=n*k,seed=seed)(N)).to(device)
    q = tf_On_QR(u,n,k)
    return q

rand_Un_QR 

rand_Un_QR(N, n, k, seed=None, device=None, qp_unif_gen=None)

Generate a batch of N random unitary matrices of size (N, n, k) using the QR decomposition.

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Number of rows.

 required
k int

Number of columns with k<=n.

 required
seed int

Random seed for reproducibility.

 None
device str

Device to store the tensor on (e.g., "cpu", "cuda").

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None

Returns:

Name Type Description
x Tensor

A batch of N random unitary matrices of size (N, n, k).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> q = rand_Un_QR(2,3,3,seed=7)
>>> q.shape 
torch.Size([2, 3, 3])
>>> q
tensor([[[-0.7005+0.0905j,  0.0129-0.0897j,  0.5259-0.4652j],
         [ 0.1364-0.4500j,  0.2766-0.6973j, -0.2259-0.4062j],
         [-0.4493+0.2795j,  0.6275-0.1878j, -0.2937+0.4525j]],

        [[ 0.5150+0.6182j, -0.1275-0.0258j, -0.5254-0.2442j],
         [ 0.1241+0.1499j, -0.6573-0.5715j,  0.3584+0.2739j],
         [ 0.1027-0.5515j, -0.0820-0.4666j, -0.2046-0.6473j]]])

>>> q = rand_Un_QR(2,3,3,seed=7,qp_unif_gen=qp.Net)
>>> q.shape 
torch.Size([2, 3, 3])
>>> q
tensor([[[ 0.1884+0.1952j,  0.2403+0.7558j,  0.1686-0.5187j],
         [-0.0775+0.3625j, -0.0203-0.5851j,  0.2945-0.6581j],
         [-0.8250-0.3293j, -0.0804+0.1475j, -0.2682-0.3329j]],

        [[-0.6121-0.3989j, -0.2673-0.1013j, -0.6161-0.0701j],
         [ 0.4590-0.0060j,  0.0882-0.5607j, -0.4561+0.5091j],
         [ 0.1663+0.4774j, -0.7245-0.2670j, -0.0730-0.3782j]]])
Source code in rmatqmcpy/rand_On_Un_Spn.py 
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
def rand_Un_QR(N, n, k, seed=None, device=None, qp_unif_gen=None):
    r"""
    Generate a batch of `N` random unitary matrices of size `(N, n, k)` using the QR decomposition. 

    Args:
        N (int): Number of samples to generate.
        n (int): Number of rows.
        k (int): Number of columns with `k<=n`.
        seed (int, optional): Random seed for reproducibility.
        device (str, optional): Device to store the tensor on (e.g., "cpu", "cuda").
        qp_unif_gen (qmcpy.DiscreteDistribution, optional): QMCPy distribution generator.

    Returns:
        x (torch.Tensor): A batch of `N` random unitary matrices of size `(N, n, k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> q = rand_Un_QR(2,3,3,seed=7)
        >>> q.shape 
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[-0.7005+0.0905j,  0.0129-0.0897j,  0.5259-0.4652j],
                 [ 0.1364-0.4500j,  0.2766-0.6973j, -0.2259-0.4062j],
                 [-0.4493+0.2795j,  0.6275-0.1878j, -0.2937+0.4525j]],
        <BLANKLINE>
                [[ 0.5150+0.6182j, -0.1275-0.0258j, -0.5254-0.2442j],
                 [ 0.1241+0.1499j, -0.6573-0.5715j,  0.3584+0.2739j],
                 [ 0.1027-0.5515j, -0.0820-0.4666j, -0.2046-0.6473j]]])

        >>> q = rand_Un_QR(2,3,3,seed=7,qp_unif_gen=qp.Net)
        >>> q.shape 
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[ 0.1884+0.1952j,  0.2403+0.7558j,  0.1686-0.5187j],
                 [-0.0775+0.3625j, -0.0203-0.5851j,  0.2945-0.6581j],
                 [-0.8250-0.3293j, -0.0804+0.1475j, -0.2682-0.3329j]],
        <BLANKLINE>
                [[-0.6121-0.3989j, -0.2673-0.1013j, -0.6161-0.0701j],
                 [ 0.4590-0.0060j,  0.0882-0.5607j, -0.4561+0.5091j],
                 [ 0.1663+0.4774j, -0.7245-0.2670j, -0.0730-0.3782j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    rng = agsutil.get_torch_rng(seed,device=device)
    u = torch.from_numpy(qp_unif_gen(dimension=2*n*k,seed=seed)(N)).to(device)
    q = tf_Un_QR(u,n,k)
    return q

rand_On_eig 

rand_On_eig(N, n, k, seed=None, device=None, qp_unif_gen=None)

Generate a batch of N random orthogonal matrices of size (N, n, k) using the eigenvalue decomposition.

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Number of rows.

 required
k int

Number of columns with k<=n.

 required
seed int

Random seed for reproducibility.

 None
device str

Device to store the tensor on (e.g., "cpu", "cuda").

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None

Returns:

Name Type Description
x Tensor

A batch of N random orthogonal matrices of size (N, n, k).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> q = rand_On_eig(2,3,3,seed=7)
>>> q.shape 
torch.Size([2, 3, 3])
>>> q
tensor([[[ 0.3748,  0.4469, -0.8123],
         [-0.1178, -0.8461, -0.5199],
         [ 0.9196, -0.2905,  0.2645]],

        [[ 0.9026, -0.2322, -0.3624],
         [-0.1140, -0.9409,  0.3190],
         [ 0.4150,  0.2467,  0.8757]]])

>>> q = rand_On_eig(2,3,3,seed=7,qp_unif_gen=qp.Net)
>>> q.shape 
torch.Size([2, 3, 3])
>>> q
tensor([[[-0.7446, -0.3508,  0.5679],
         [-0.5217,  0.8366, -0.1671],
         [-0.4165, -0.4207, -0.8060]],

        [[-0.0434, -0.8815,  0.4702],
         [-0.2524,  0.4650,  0.8486],
         [ 0.9666,  0.0818,  0.2427]]])
Source code in rmatqmcpy/rand_On_Un_Spn.py 
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
def rand_On_eig(N, n, k, seed=None, device=None, qp_unif_gen=None):
    r"""
    Generate a batch of `N` random orthogonal matrices of size `(N, n, k)` using the eigenvalue decomposition.

    Args:
        N (int): Number of samples to generate.
        n (int): Number of rows.
        k (int): Number of columns with `k<=n`.
        seed (int, optional): Random seed for reproducibility.
        device (str, optional): Device to store the tensor on (e.g., "cpu", "cuda").
        qp_unif_gen (qmcpy.DiscreteDistribution, optional): QMCPy distribution generator.

    Returns:
        x (torch.Tensor): A batch of `N` random orthogonal matrices of size `(N, n, k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> q = rand_On_eig(2,3,3,seed=7)
        >>> q.shape 
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[ 0.3748,  0.4469, -0.8123],
                 [-0.1178, -0.8461, -0.5199],
                 [ 0.9196, -0.2905,  0.2645]],
        <BLANKLINE>
                [[ 0.9026, -0.2322, -0.3624],
                 [-0.1140, -0.9409,  0.3190],
                 [ 0.4150,  0.2467,  0.8757]]])

        >>> q = rand_On_eig(2,3,3,seed=7,qp_unif_gen=qp.Net)
        >>> q.shape 
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[-0.7446, -0.3508,  0.5679],
                 [-0.5217,  0.8366, -0.1671],
                 [-0.4165, -0.4207, -0.8060]],
        <BLANKLINE>
                [[-0.0434, -0.8815,  0.4702],
                 [-0.2524,  0.4650,  0.8486],
                 [ 0.9666,  0.0818,  0.2427]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    u = torch.from_numpy(qp_unif_gen(dimension=n*(n-1)//2+2*n,seed=seed)(N)).to(device)
    q = tf_On_eig(u,n,k)
    return q

rand_Un_eig 

rand_Un_eig(N, n, k, seed=None, device=None, qp_unif_gen=None)

Generate a batch of N random unitary matrices of size (N, n, k) using the eigenvalue decomposition.

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Number of rows.

 required
k int

Number of columns with k<=n.

 required
seed int

Random seed for reproducibility.

 None
device str

Device to store the tensor on (e.g., "cpu", "cuda").

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None

Returns:

Name Type Description
x Tensor

A batch of N random unitary matrices of size (N, n, k).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> q = rand_Un_eig(2,3,3,seed=7)
>>> q.shape 
torch.Size([2, 3, 3])
>>> q
tensor([[[-0.4009-0.1134j,  0.5152+0.3147j, -0.4627+0.4978j],
         [-0.0383-0.0098j, -0.5419+0.5188j,  0.3717+0.5453j],
         [-0.5107-0.7510j, -0.0557-0.2638j,  0.3130-0.0670j]],

        [[ 0.1113+0.1290j,  0.0108-0.7796j, -0.3388+0.4983j],
         [ 0.0135-0.9445j,  0.1002+0.0782j, -0.2065+0.2213j],
         [ 0.1714+0.2221j,  0.3955+0.4685j,  0.1415+0.7248j]]])

>>> q = rand_Un_eig(2,3,3,seed=7,qp_unif_gen=qp.Net)
>>> q.shape 
torch.Size([2, 3, 3])
>>> q
tensor([[[ 0.5988-0.2838j, -0.3694+0.2998j, -0.5360+0.2174j],
         [ 0.0116-0.3634j,  0.6324-0.2327j, -0.5162-0.3838j],
         [-0.6395+0.1405j, -0.2209+0.5204j, -0.4418-0.2376j]],

        [[-0.6903+0.0581j, -0.1665+0.1366j, -0.4519-0.5192j],
         [-0.6830+0.0339j,  0.0105-0.3853j,  0.3131+0.5346j],
         [-0.2250-0.0432j,  0.2931+0.8480j,  0.3613+0.1091j]]])
Source code in rmatqmcpy/rand_On_Un_Spn.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
def rand_Un_eig(N, n, k, seed=None, device=None, qp_unif_gen=None):
    r"""
    Generate a batch of `N` random unitary matrices of size `(N, n, k)` using the eigenvalue decomposition.

    Args:
        N (int): Number of samples to generate.
        n (int): Number of rows.
        k (int): Number of columns with `k<=n`.
        seed (int, optional): Random seed for reproducibility.
        device (str, optional): Device to store the tensor on (e.g., "cpu", "cuda").
        qp_unif_gen (qmcpy.DiscreteDistribution, optional): QMCPy distribution generator.

    Returns:
        x (torch.Tensor): A batch of `N` random unitary matrices of size `(N, n, k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> q = rand_Un_eig(2,3,3,seed=7)
        >>> q.shape 
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[-0.4009-0.1134j,  0.5152+0.3147j, -0.4627+0.4978j],
                 [-0.0383-0.0098j, -0.5419+0.5188j,  0.3717+0.5453j],
                 [-0.5107-0.7510j, -0.0557-0.2638j,  0.3130-0.0670j]],
        <BLANKLINE>
                [[ 0.1113+0.1290j,  0.0108-0.7796j, -0.3388+0.4983j],
                 [ 0.0135-0.9445j,  0.1002+0.0782j, -0.2065+0.2213j],
                 [ 0.1714+0.2221j,  0.3955+0.4685j,  0.1415+0.7248j]]])

        >>> q = rand_Un_eig(2,3,3,seed=7,qp_unif_gen=qp.Net)
        >>> q.shape 
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[ 0.5988-0.2838j, -0.3694+0.2998j, -0.5360+0.2174j],
                 [ 0.0116-0.3634j,  0.6324-0.2327j, -0.5162-0.3838j],
                 [-0.6395+0.1405j, -0.2209+0.5204j, -0.4418-0.2376j]],
        <BLANKLINE>
                [[-0.6903+0.0581j, -0.1665+0.1366j, -0.4519-0.5192j],
                 [-0.6830+0.0339j,  0.0105-0.3853j,  0.3131+0.5346j],
                 [-0.2250-0.0432j,  0.2931+0.8480j,  0.3613+0.1091j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    u = torch.from_numpy(qp_unif_gen(dimension=n**2+n,seed=seed)(N)).to(device)
    q = tf_Un_eig(u,n,k)
    return q

rand_Spn_SVD 

rand_Spn_SVD(N, n, k, seed=None, device=None, qp_unif_gen=None)

Generate a batch of N random symplectic unitary matrices of size (N, 2n, 2k) using the SVD decomposition.

Parameters:

Name Type Description Default
N int

Number of samples to generate.

 required
n int

Half the number of rows.

 required
k int

Half the number of columns with k<=n.

 required
seed int

Random seed for reproducibility.

 None
device str

Device to store the tensor on (e.g., "cpu", "cuda").

 None
qp_unif_gen DiscreteDistribution

QMCPy distribution generator.

 None

Returns:

Name Type Description
x Tensor

A batch of N random symplectic unitary matrices of size (N, 2n, 2k).

Examples:

>>> torch.set_default_dtype(torch.float64) 

>>> q = rand_Spn_SVD(2,3,3,seed=7)
>>> q.shape 
torch.Size([2, 6, 6])
>>> q
tensor([[[-0.5379+0.0018j,  0.1665-0.3441j,  0.1465-0.4183j, -0.0826+0.2936j,
           0.0846-0.3585j, -0.3592-0.1015j],
         [-0.2501+0.3772j, -0.4710+0.0421j,  0.1512+0.1012j,  0.1678-0.1855j,
           0.3810+0.2695j, -0.1729-0.4776j],
         [ 0.0476+0.3640j,  0.4540-0.1432j, -0.1454-0.3261j,  0.1149-0.4497j,
           0.1691-0.1465j,  0.4496-0.2084j],
         [ 0.0826+0.2936j, -0.0846-0.3585j,  0.3592-0.1015j, -0.5379-0.0018j,
           0.1665+0.3441j,  0.1465+0.4183j],
         [-0.1678-0.1855j, -0.3810+0.2695j,  0.1729-0.4776j, -0.2501-0.3772j,
          -0.4710-0.0421j,  0.1512-0.1012j],
         [-0.1149-0.4497j, -0.1691-0.1465j, -0.4496-0.2084j,  0.0476-0.3640j,
           0.4540+0.1432j, -0.1454+0.3261j]],

        [[-0.0798-0.5452j, -0.1923+0.2617j,  0.0815+0.2038j, -0.0359-0.4116j,
           0.2102-0.3235j,  0.4459+0.1562j],
         [-0.3699+0.0262j, -0.2731-0.5066j, -0.1090-0.0730j, -0.2111-0.4440j,
          -0.1848-0.1637j, -0.3600+0.2859j],
         [ 0.1009+0.1924j,  0.1127-0.3529j,  0.5547+0.1064j,  0.2102+0.2346j,
          -0.2066-0.4166j,  0.2735+0.3259j],
         [ 0.0359-0.4116j, -0.2102-0.3235j, -0.4459+0.1562j, -0.0798+0.5452j,
          -0.1923-0.2617j,  0.0815-0.2038j],
         [ 0.2111-0.4440j,  0.1848-0.1637j,  0.3600+0.2859j, -0.3699-0.0262j,
          -0.2731+0.5066j, -0.1090+0.0730j],
         [-0.2102+0.2346j,  0.2066-0.4166j, -0.2735+0.3259j,  0.1009-0.1924j,
           0.1127+0.3529j,  0.5547-0.1064j]]])

>>> q = rand_Spn_SVD(2,3,3,seed=7,qp_unif_gen=qp.Net)
>>> q.shape 
torch.Size([2, 6, 6])
>>> q
tensor([[[ 0.0024+0.0842j,  0.6149-0.0634j, -0.3828-0.0823j, -0.2519-0.3091j,
           0.3628+0.2602j, -0.1177+0.2920j],
         [ 0.3461-0.6176j,  0.0761+0.1610j, -0.5182+0.0585j,  0.2130+0.0398j,
          -0.1930-0.0418j,  0.3297+0.0230j],
         [ 0.0650-0.1395j,  0.2626-0.2246j,  0.3886-0.2880j,  0.4926+0.1394j,
           0.4424+0.1808j,  0.3377-0.1361j],
         [ 0.2519-0.3091j, -0.3628+0.2602j,  0.1177+0.2920j,  0.0024-0.0842j,
           0.6149+0.0634j, -0.3828+0.0823j],
         [-0.2130+0.0398j,  0.1930-0.0418j, -0.3297+0.0230j,  0.3461+0.6176j,
           0.0761-0.1610j, -0.5182-0.0585j],
         [-0.4926+0.1394j, -0.4424+0.1808j, -0.3377-0.1361j,  0.0650+0.1395j,
           0.2626+0.2246j,  0.3886+0.2880j]],

        [[ 0.2988+0.4254j, -0.3412+0.0841j,  0.2553+0.2561j,  0.0461+0.1016j,
          -0.1127-0.1724j,  0.5693-0.3106j],
         [-0.1268+0.1878j,  0.2846+0.1790j,  0.2376+0.0494j, -0.3145-0.5206j,
           0.0981-0.6136j, -0.1342-0.0512j],
         [ 0.4454-0.0571j, -0.0687+0.2563j, -0.5844+0.0523j, -0.0528+0.3026j,
           0.2235-0.4632j, -0.1256+0.0951j],
         [-0.0461+0.1016j,  0.1127-0.1724j, -0.5693-0.3106j,  0.2988-0.4254j,
          -0.3412-0.0841j,  0.2553-0.2561j],
         [ 0.3145-0.5206j, -0.0981-0.6136j,  0.1342-0.0512j, -0.1268-0.1878j,
           0.2846-0.1790j,  0.2376-0.0494j],
         [ 0.0528+0.3026j, -0.2235-0.4632j,  0.1256+0.0951j,  0.4454+0.0571j,
          -0.0687-0.2563j, -0.5844-0.0523j]]])
Source code in rmatqmcpy/rand_On_Un_Spn.py 
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
def rand_Spn_SVD(N, n, k, seed=None, device=None, qp_unif_gen=None):
    r"""
    Generate a batch of `N` random symplectic unitary matrices of size `(N, 2n, 2k)` using the SVD decomposition. 

    Args:
        N (int): Number of samples to generate.
        n (int): Half the number of rows.
        k (int): Half the number of columns with `k<=n`.
        seed (int, optional): Random seed for reproducibility.
        device (str, optional): Device to store the tensor on (e.g., "cpu", "cuda").
        qp_unif_gen (qmcpy.DiscreteDistribution, optional): QMCPy distribution generator.

    Returns:
        x (torch.Tensor): A batch of `N` random symplectic unitary matrices of size `(N, 2n, 2k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 

        >>> q = rand_Spn_SVD(2,3,3,seed=7)
        >>> q.shape 
        torch.Size([2, 6, 6])
        >>> q
        tensor([[[-0.5379+0.0018j,  0.1665-0.3441j,  0.1465-0.4183j, -0.0826+0.2936j,
                   0.0846-0.3585j, -0.3592-0.1015j],
                 [-0.2501+0.3772j, -0.4710+0.0421j,  0.1512+0.1012j,  0.1678-0.1855j,
                   0.3810+0.2695j, -0.1729-0.4776j],
                 [ 0.0476+0.3640j,  0.4540-0.1432j, -0.1454-0.3261j,  0.1149-0.4497j,
                   0.1691-0.1465j,  0.4496-0.2084j],
                 [ 0.0826+0.2936j, -0.0846-0.3585j,  0.3592-0.1015j, -0.5379-0.0018j,
                   0.1665+0.3441j,  0.1465+0.4183j],
                 [-0.1678-0.1855j, -0.3810+0.2695j,  0.1729-0.4776j, -0.2501-0.3772j,
                  -0.4710-0.0421j,  0.1512-0.1012j],
                 [-0.1149-0.4497j, -0.1691-0.1465j, -0.4496-0.2084j,  0.0476-0.3640j,
                   0.4540+0.1432j, -0.1454+0.3261j]],
        <BLANKLINE>
                [[-0.0798-0.5452j, -0.1923+0.2617j,  0.0815+0.2038j, -0.0359-0.4116j,
                   0.2102-0.3235j,  0.4459+0.1562j],
                 [-0.3699+0.0262j, -0.2731-0.5066j, -0.1090-0.0730j, -0.2111-0.4440j,
                  -0.1848-0.1637j, -0.3600+0.2859j],
                 [ 0.1009+0.1924j,  0.1127-0.3529j,  0.5547+0.1064j,  0.2102+0.2346j,
                  -0.2066-0.4166j,  0.2735+0.3259j],
                 [ 0.0359-0.4116j, -0.2102-0.3235j, -0.4459+0.1562j, -0.0798+0.5452j,
                  -0.1923-0.2617j,  0.0815-0.2038j],
                 [ 0.2111-0.4440j,  0.1848-0.1637j,  0.3600+0.2859j, -0.3699-0.0262j,
                  -0.2731+0.5066j, -0.1090+0.0730j],
                 [-0.2102+0.2346j,  0.2066-0.4166j, -0.2735+0.3259j,  0.1009-0.1924j,
                   0.1127+0.3529j,  0.5547-0.1064j]]])

        >>> q = rand_Spn_SVD(2,3,3,seed=7,qp_unif_gen=qp.Net)
        >>> q.shape 
        torch.Size([2, 6, 6])
        >>> q
        tensor([[[ 0.0024+0.0842j,  0.6149-0.0634j, -0.3828-0.0823j, -0.2519-0.3091j,
                   0.3628+0.2602j, -0.1177+0.2920j],
                 [ 0.3461-0.6176j,  0.0761+0.1610j, -0.5182+0.0585j,  0.2130+0.0398j,
                  -0.1930-0.0418j,  0.3297+0.0230j],
                 [ 0.0650-0.1395j,  0.2626-0.2246j,  0.3886-0.2880j,  0.4926+0.1394j,
                   0.4424+0.1808j,  0.3377-0.1361j],
                 [ 0.2519-0.3091j, -0.3628+0.2602j,  0.1177+0.2920j,  0.0024-0.0842j,
                   0.6149+0.0634j, -0.3828+0.0823j],
                 [-0.2130+0.0398j,  0.1930-0.0418j, -0.3297+0.0230j,  0.3461+0.6176j,
                   0.0761-0.1610j, -0.5182-0.0585j],
                 [-0.4926+0.1394j, -0.4424+0.1808j, -0.3377-0.1361j,  0.0650+0.1395j,
                   0.2626+0.2246j,  0.3886+0.2880j]],
        <BLANKLINE>
                [[ 0.2988+0.4254j, -0.3412+0.0841j,  0.2553+0.2561j,  0.0461+0.1016j,
                  -0.1127-0.1724j,  0.5693-0.3106j],
                 [-0.1268+0.1878j,  0.2846+0.1790j,  0.2376+0.0494j, -0.3145-0.5206j,
                   0.0981-0.6136j, -0.1342-0.0512j],
                 [ 0.4454-0.0571j, -0.0687+0.2563j, -0.5844+0.0523j, -0.0528+0.3026j,
                   0.2235-0.4632j, -0.1256+0.0951j],
                 [-0.0461+0.1016j,  0.1127-0.1724j, -0.5693-0.3106j,  0.2988-0.4254j,
                  -0.3412-0.0841j,  0.2553-0.2561j],
                 [ 0.3145-0.5206j, -0.0981-0.6136j,  0.1342-0.0512j, -0.1268-0.1878j,
                   0.2846-0.1790j,  0.2376-0.0494j],
                 [ 0.0528+0.3026j, -0.2235-0.4632j,  0.1256+0.0951j,  0.4454+0.0571j,
                  -0.0687-0.2563j, -0.5844-0.0523j]]])
    """
    if qp_unif_gen is None: qp_unif_gen = qp.IIDStdUniform
    u = torch.from_numpy(qp_unif_gen(dimension=4*n**2,seed=seed)(N)).to(device)
    q = tf_Spn_SVD(u,n,k)
    return q

Transforms from uniforms to \(\mathrm{O}(n)\), \(\mathrm{U}(n)\) and \(\mathrm{Sp}(n)\)

tf_On_QR 

tf_On_QR(u, n, k)

Transform uniform samples to random orthogonal matrices using QR decomposition.

Parameters:

Name Type Description Default
u Tensor

Uniformly distributed samples in [0, 1] of shape (..., n*k).

 required
n int

Number of rows.

 required
k int

Number of columns with k<=n.

 required

Returns:

Name Type Description
x Tensor

A batch of random orthogonal matrices of shape (..., n, k).

Examples:

>>> torch.set_default_dtype(torch.float64) 
>>> rng = torch.Generator().manual_seed(7)

>>> delta = torch.arange(1,4,dtype=float)
>>> delta = delta/torch.linalg.norm(delta)
>>> n = delta.size(-1)
>>> q = tf_On_QR(torch.rand((2,n**2),generator=rng),n,n)
>>> q.shape
torch.Size([2, 3, 3])
>>> q
tensor([[[-0.6554,  0.0144,  0.7552],
         [ 0.4522,  0.8083,  0.3771],
         [-0.6049,  0.5886, -0.5363]],

        [[ 0.4473, -0.3480, -0.8239],
         [-0.6125, -0.7904,  0.0013],
         [ 0.6517, -0.5041,  0.5667]]])
>>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q),torch.eye(n))
True

>>> k = 2
>>> q = tf_On_QR(torch.rand((4,n*k),generator=rng),n,k)
>>> q.shape
torch.Size([4, 3, 2])
>>> q
tensor([[[-0.6351,  0.7368],
         [ 0.7422,  0.6652],
         [ 0.2139, -0.1206]],

        [[ 0.5383,  0.2390],
         [ 0.7023,  0.4028],
         [-0.4658,  0.8835]],

        [[ 0.1407, -0.9139],
         [ 0.6684, -0.1959],
         [-0.7303, -0.3554]],

        [[-0.7980, -0.4883],
         [ 0.5758, -0.4442],
         [-0.1782,  0.7512]]])
>>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q),torch.eye(k))
True
Source code in rmatqmcpy/tf_On_Un_Spn.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def tf_On_QR(u, n, k):
    r"""
    Transform uniform samples to random orthogonal matrices using QR decomposition.

    Args:
        u (torch.Tensor): Uniformly distributed samples in [0, 1] of shape `(..., n*k)`.
        n (int): Number of rows.
        k (int): Number of columns with `k<=n`.

    Returns:
        x (torch.Tensor): A batch of random orthogonal matrices of shape `(..., n, k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 
        >>> rng = torch.Generator().manual_seed(7)

        >>> delta = torch.arange(1,4,dtype=float)
        >>> delta = delta/torch.linalg.norm(delta)
        >>> n = delta.size(-1)
        >>> q = tf_On_QR(torch.rand((2,n**2),generator=rng),n,n)
        >>> q.shape
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[-0.6554,  0.0144,  0.7552],
                 [ 0.4522,  0.8083,  0.3771],
                 [-0.6049,  0.5886, -0.5363]],
        <BLANKLINE>
                [[ 0.4473, -0.3480, -0.8239],
                 [-0.6125, -0.7904,  0.0013],
                 [ 0.6517, -0.5041,  0.5667]]])
        >>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q),torch.eye(n))
        True

        >>> k = 2
        >>> q = tf_On_QR(torch.rand((4,n*k),generator=rng),n,k)
        >>> q.shape
        torch.Size([4, 3, 2])
        >>> q
        tensor([[[-0.6351,  0.7368],
                 [ 0.7422,  0.6652],
                 [ 0.2139, -0.1206]],
        <BLANKLINE>
                [[ 0.5383,  0.2390],
                 [ 0.7023,  0.4028],
                 [-0.4658,  0.8835]],
        <BLANKLINE>
                [[ 0.1407, -0.9139],
                 [ 0.6684, -0.1959],
                 [-0.7303, -0.3554]],
        <BLANKLINE>
                [[-0.7980, -0.4883],
                 [ 0.5758, -0.4442],
                 [-0.1782,  0.7512]]])
        >>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q),torch.eye(k))
        True
    """ 
    assert k<=n
    assert u.size(-1)==(n*k)
    assert (0<=u).all()
    assert (u<=1).all()
    u = u.reshape((*u.shape[:-1],n,k))
    v = agsutil.icdf_std_normal(u)
    q,r = torch.linalg.qr(v,mode='reduced')
    d = torch.diagonal(r,dim1=-2,dim2=-1).sign()
    q = q*d[...,None,:]
    # assert torch.allclose(torch.einsum("...ij,...jk->...ik",q,r*d[...,:,None]),v)
    return q

tf_Un_QR 

tf_Un_QR(u, n, k)

Transform uniform samples to random unitary matrices using QR decomposition.

Parameters:

Name Type Description Default
u Tensor

Uniformly distributed samples in [0, 1] of shape (..., 2*n*k).

 required
n int

Number of rows.

 required
k int

Number of columns with k<=n.

 required

Returns:

Name Type Description
x Tensor

A batch of random unitary matrices of shape (..., n, k).

Examples:

>>> torch.set_default_dtype(torch.float64) 
>>> rng = torch.Generator().manual_seed(7)

>>> delta = torch.arange(1,4,dtype=float)
>>> delta = delta/torch.linalg.norm(delta)
>>> n = delta.size(-1)
>>> q = tf_Un_QR(torch.rand((2,2*n**2),generator=rng),n,n)
>>> q.shape
torch.Size([2, 3, 3])
>>> q
tensor([[[-0.2586-0.2661j,  0.8447-0.1177j,  0.3613+0.0667j],
         [-0.2387+0.0676j,  0.2890+0.2913j, -0.5804-0.6582j],
         [-0.4265-0.7867j, -0.2891-0.1438j, -0.2795+0.1297j]],

        [[-0.3620+0.7348j,  0.3715+0.0140j,  0.3549-0.2545j],
         [ 0.3698+0.1818j,  0.4820+0.5669j, -0.1691+0.4979j],
         [ 0.0229-0.3983j,  0.3680+0.4155j, -0.0743-0.7261j]]])
>>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(n,dtype=q.dtype))
True

>>> k = 2
>>> q = tf_Un_QR(torch.rand((4,2*n*k),generator=rng),n,k)
>>> q.shape
torch.Size([4, 3, 2])
>>> q
tensor([[[-0.3500-0.6164j, -0.2414-0.1013j],
         [-0.0782+0.5104j,  0.4290+0.1885j],
         [-0.3047-0.3715j,  0.8351-0.1207j]],

        [[ 0.0614+0.3318j, -0.8324-0.2894j],
         [-0.3029+0.5080j,  0.3375+0.3046j],
         [ 0.6301+0.3731j,  0.1075+0.0721j]],

        [[ 0.7912+0.2576j,  0.4185+0.2187j],
         [ 0.3229+0.0629j, -0.7973-0.0245j],
         [ 0.0592-0.4425j, -0.2767+0.2533j]],

        [[-0.2778-0.3864j, -0.2751+0.8253j],
         [-0.2274+0.0779j,  0.0250-0.1292j],
         [ 0.7963+0.2856j,  0.1643+0.4459j]]])
>>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(k,dtype=q.dtype))
True
Source code in rmatqmcpy/tf_On_Un_Spn.py 
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
def tf_Un_QR(u, n, k):
    r"""
    Transform uniform samples to random unitary matrices using QR decomposition.

    Args:
        u (torch.Tensor): Uniformly distributed samples in [0, 1] of shape `(..., 2*n*k)`.
        n (int): Number of rows.
        k (int): Number of columns with `k<=n`.

    Returns:
        x (torch.Tensor): A batch of random unitary matrices of shape `(..., n, k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 
        >>> rng = torch.Generator().manual_seed(7)

        >>> delta = torch.arange(1,4,dtype=float)
        >>> delta = delta/torch.linalg.norm(delta)
        >>> n = delta.size(-1)
        >>> q = tf_Un_QR(torch.rand((2,2*n**2),generator=rng),n,n)
        >>> q.shape
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[-0.2586-0.2661j,  0.8447-0.1177j,  0.3613+0.0667j],
                 [-0.2387+0.0676j,  0.2890+0.2913j, -0.5804-0.6582j],
                 [-0.4265-0.7867j, -0.2891-0.1438j, -0.2795+0.1297j]],
        <BLANKLINE>
                [[-0.3620+0.7348j,  0.3715+0.0140j,  0.3549-0.2545j],
                 [ 0.3698+0.1818j,  0.4820+0.5669j, -0.1691+0.4979j],
                 [ 0.0229-0.3983j,  0.3680+0.4155j, -0.0743-0.7261j]]])
        >>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(n,dtype=q.dtype))
        True

        >>> k = 2
        >>> q = tf_Un_QR(torch.rand((4,2*n*k),generator=rng),n,k)
        >>> q.shape
        torch.Size([4, 3, 2])
        >>> q
        tensor([[[-0.3500-0.6164j, -0.2414-0.1013j],
                 [-0.0782+0.5104j,  0.4290+0.1885j],
                 [-0.3047-0.3715j,  0.8351-0.1207j]],
        <BLANKLINE>
                [[ 0.0614+0.3318j, -0.8324-0.2894j],
                 [-0.3029+0.5080j,  0.3375+0.3046j],
                 [ 0.6301+0.3731j,  0.1075+0.0721j]],
        <BLANKLINE>
                [[ 0.7912+0.2576j,  0.4185+0.2187j],
                 [ 0.3229+0.0629j, -0.7973-0.0245j],
                 [ 0.0592-0.4425j, -0.2767+0.2533j]],
        <BLANKLINE>
                [[-0.2778-0.3864j, -0.2751+0.8253j],
                 [-0.2274+0.0779j,  0.0250-0.1292j],
                 [ 0.7963+0.2856j,  0.1643+0.4459j]]])
        >>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(k,dtype=q.dtype))
        True
    """ 
    assert k<=n
    assert u.size(-1)==(2*n*k)
    assert (0<=u).all()
    assert (u<=1).all()
    u = u.reshape((*u.shape[:-1],n,k,2))
    v = agsutil.icdf_std_normal(u)
    v_complex = torch.complex(v[...,0],v[...,1])
    q,r = torch.linalg.qr(v_complex,mode='reduced')
    d = torch.diagonal(r,dim1=-2,dim2=-1)
    ph = torch.sgn(d)
    q = q*ph[...,None,:]
    # assert torch.allclose(torch.einsum("...ij,...jk->...ik",q,r*ph[..., :, None].conj()),v_complex)
    return q

tf_On_eig 

tf_On_eig(u, n, k)

Transform uniform samples to random orthogonal matrices using eigenvalue decomposition.

Parameters:

Name Type Description Default
u Tensor

Uniformly distributed samples in [0, 1] of shape (..., n*(n+1)//2 + 2*n).

 required
n int

Number of rows.

 required
k int

Number of columns with k<=n.

 required

Returns:

Name Type Description
x Tensor

A batch of random orthogonal matrices of shape (..., n, k).

Examples:

>>> torch.set_default_dtype(torch.float64) 
>>> rng = torch.Generator().manual_seed(7)

>>> delta = torch.arange(1,4,dtype=float)
>>> delta = delta/torch.linalg.norm(delta)
>>> n = delta.size(-1)
>>> q = tf_On_eig(torch.rand((2,2*n+n*(n-1)//2),generator=rng),n,n)
>>> q.shape
torch.Size([2, 3, 3])
>>> q
tensor([[[-0.8410,  0.5062, -0.1912],
         [-0.4034, -0.3511,  0.8450],
         [ 0.3606,  0.7877,  0.4995]],

        [[ 0.5232, -0.5182,  0.6765],
         [-0.5917, -0.7922, -0.1492],
         [ 0.6133, -0.3223, -0.7211]]])
>>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q),torch.eye(n))
True

>>> k = 2
>>> q = tf_On_eig(torch.rand((4,2*n+n*(n-1)//2),generator=rng),n,k)
>>> q.shape
torch.Size([4, 3, 2])
>>> q
tensor([[[ 0.0680,  0.7201],
         [-0.5824, -0.5334],
         [ 0.8101, -0.4439]],

        [[ 0.1218, -0.9108],
         [-0.8979,  0.0683],
         [-0.4230, -0.4071]],

        [[-0.7189,  0.6946],
         [-0.6655, -0.6771],
         [-0.2007, -0.2431]],

        [[ 0.1882, -0.9249],
         [-0.2634, -0.3716],
         [ 0.9461,  0.0805]]])
>>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q),torch.eye(k))
True
Source code in rmatqmcpy/tf_On_Un_Spn.py 
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def tf_On_eig(u, n, k):
    r"""
    Transform uniform samples to random orthogonal matrices  using eigenvalue decomposition.

    Args:
        u (torch.Tensor): Uniformly distributed samples in [0, 1] of shape `(..., n*(n+1)//2 + 2*n)`.
        n (int): Number of rows.
        k (int): Number of columns with `k<=n`.

    Returns:
        x (torch.Tensor): A batch of random orthogonal matrices of shape `(..., n, k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 
        >>> rng = torch.Generator().manual_seed(7)

        >>> delta = torch.arange(1,4,dtype=float)
        >>> delta = delta/torch.linalg.norm(delta)
        >>> n = delta.size(-1)
        >>> q = tf_On_eig(torch.rand((2,2*n+n*(n-1)//2),generator=rng),n,n)
        >>> q.shape
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[-0.8410,  0.5062, -0.1912],
                 [-0.4034, -0.3511,  0.8450],
                 [ 0.3606,  0.7877,  0.4995]],
        <BLANKLINE>
                [[ 0.5232, -0.5182,  0.6765],
                 [-0.5917, -0.7922, -0.1492],
                 [ 0.6133, -0.3223, -0.7211]]])
        >>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q),torch.eye(n))
        True

        >>> k = 2
        >>> q = tf_On_eig(torch.rand((4,2*n+n*(n-1)//2),generator=rng),n,k)
        >>> q.shape
        torch.Size([4, 3, 2])
        >>> q
        tensor([[[ 0.0680,  0.7201],
                 [-0.5824, -0.5334],
                 [ 0.8101, -0.4439]],
        <BLANKLINE>
                [[ 0.1218, -0.9108],
                 [-0.8979,  0.0683],
                 [-0.4230, -0.4071]],
        <BLANKLINE>
                [[-0.7189,  0.6946],
                 [-0.6655, -0.6771],
                 [-0.2007, -0.2431]],
        <BLANKLINE>
                [[ 0.1882, -0.9249],
                 [-0.2634, -0.3716],
                 [ 0.9461,  0.0805]]])
        >>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q),torch.eye(k))
        True
    """
    assert k<=n
    num_tril = n*(n-1)//2
    assert u.size(-1)==(2*n+num_tril)
    assert (0<=u).all()
    assert (u<=1).all()
    signs = torch.where(u[...,:n]>1/2,1.,-1.).to(u.device)
    alpha = agsutil.icdf_std_normal(u[...,n:2*n])*np.sqrt(2)
    beta = agsutil.icdf_std_normal(u[...,2*n:2*n+num_tril])
    il0,il1 = torch.tril_indices(n,n,offset=-1,device=u.device)
    v = torch.zeros((*u.shape[:-1],n,n),dtype=u.dtype,device=u.device)
    diag_indices = torch.arange(n,device=u.device)
    v[...,diag_indices,diag_indices] = alpha
    v[...,il0,il1] = beta
    v += v.tril(-1).transpose(dim0=-2,dim1=-1)
    assert torch.allclose(v[...,torch.arange(n,device=u.device),torch.arange(n,device=u.device)],alpha)
    gamma,q = torch.linalg.eigh(v)
    q = q*signs[...,None,:]
    # assert torch.allclose(torch.einsum("...ij,...j,...kj->...ik",q,gamma,q),v)
    return q[...,:k]

tf_Un_eig 

tf_Un_eig(u, n, k)

Transform uniform samples to random unitary matrices using eigenvalue decomposition.

Parameters:

Name Type Description Default
u Tensor

Uniformly distributed samples in [0, 1] of shape (..., n**2 + 2*n).

 required
n int

Number of rows.

 required
k int

Number of columns with k<=n.

 required

Returns:

Name Type Description
x Tensor

A batch of random unitary matrices of shape (..., n, k).

Examples:

>>> torch.set_default_dtype(torch.float64) 
>>> rng = torch.Generator().manual_seed(7)

>>> delta = torch.arange(1,4,dtype=float)
>>> delta = delta/torch.linalg.norm(delta)
>>> n = delta.size(-1)
>>> q = tf_Un_eig(torch.rand((2,n**2+n),generator=rng),n,n)
>>> q.shape
torch.Size([2, 3, 3])
>>> q
tensor([[[ 0.1184-0.6343j, -0.1075+0.7167j,  0.1566-0.1844j],
         [-0.1909-0.4149j, -0.0457-0.5509j,  0.6462-0.2611j],
         [ 0.6050+0.0953j,  0.3971-0.1076j, -0.0938-0.6685j]],

        [[-0.1771-0.3102j,  0.6653+0.1603j, -0.4573-0.4415j],
         [-0.7242-0.5052j, -0.4297-0.0679j,  0.1084-0.1391j],
         [ 0.2435+0.1829j, -0.4473-0.3773j, -0.3695-0.6544j]]])
>>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(n,dtype=q.dtype))
True

>>> k = 2
>>> q = tf_Un_eig(torch.rand((4,n**2+n),generator=rng),n,k)
>>> q.shape
torch.Size([4, 3, 2])
>>> q
tensor([[[-0.2158+0.1676j, -0.0992-0.8279j],
         [ 0.5745+0.6733j,  0.4101-0.0443j],
         [-0.3384+0.1658j,  0.0866-0.3565j]],

        [[ 0.2821+0.7138j, -0.4911-0.1984j],
         [-0.4658+0.4386j, -0.0056+0.6666j],
         [ 0.0355-0.0141j, -0.2300+0.4714j]],

        [[-0.4660-0.2014j, -0.2798+0.6882j],
         [-0.4888+0.5990j, -0.2285-0.4499j],
         [ 0.3796+0.0216j,  0.4329+0.0782j]],

        [[ 0.2325-0.7290j,  0.1024+0.0821j],
         [ 0.4156-0.3091j,  0.0759+0.4502j],
         [ 0.3812+0.0299j,  0.4367-0.7640j]]])
>>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(k,dtype=q.dtype))
True
Source code in rmatqmcpy/tf_On_Un_Spn.py 
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
def tf_Un_eig(u, n, k):
    r"""
    Transform uniform samples to random unitary matrices using eigenvalue decomposition.

    Args:
        u (torch.Tensor): Uniformly distributed samples in [0, 1] of shape `(..., n**2 + 2*n)`.
        n (int): Number of rows.
        k (int): Number of columns with `k<=n`.

    Returns:
        x (torch.Tensor): A batch of random unitary matrices of shape `(..., n, k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 
        >>> rng = torch.Generator().manual_seed(7)

        >>> delta = torch.arange(1,4,dtype=float)
        >>> delta = delta/torch.linalg.norm(delta)
        >>> n = delta.size(-1)
        >>> q = tf_Un_eig(torch.rand((2,n**2+n),generator=rng),n,n)
        >>> q.shape
        torch.Size([2, 3, 3])
        >>> q
        tensor([[[ 0.1184-0.6343j, -0.1075+0.7167j,  0.1566-0.1844j],
                 [-0.1909-0.4149j, -0.0457-0.5509j,  0.6462-0.2611j],
                 [ 0.6050+0.0953j,  0.3971-0.1076j, -0.0938-0.6685j]],
        <BLANKLINE>
                [[-0.1771-0.3102j,  0.6653+0.1603j, -0.4573-0.4415j],
                 [-0.7242-0.5052j, -0.4297-0.0679j,  0.1084-0.1391j],
                 [ 0.2435+0.1829j, -0.4473-0.3773j, -0.3695-0.6544j]]])
        >>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(n,dtype=q.dtype))
        True

        >>> k = 2
        >>> q = tf_Un_eig(torch.rand((4,n**2+n),generator=rng),n,k)
        >>> q.shape
        torch.Size([4, 3, 2])
        >>> q
        tensor([[[-0.2158+0.1676j, -0.0992-0.8279j],
                 [ 0.5745+0.6733j,  0.4101-0.0443j],
                 [-0.3384+0.1658j,  0.0866-0.3565j]],
        <BLANKLINE>
                [[ 0.2821+0.7138j, -0.4911-0.1984j],
                 [-0.4658+0.4386j, -0.0056+0.6666j],
                 [ 0.0355-0.0141j, -0.2300+0.4714j]],
        <BLANKLINE>
                [[-0.4660-0.2014j, -0.2798+0.6882j],
                 [-0.4888+0.5990j, -0.2285-0.4499j],
                 [ 0.3796+0.0216j,  0.4329+0.0782j]],
        <BLANKLINE>
                [[ 0.2325-0.7290j,  0.1024+0.0821j],
                 [ 0.4156-0.3091j,  0.0759+0.4502j],
                 [ 0.3812+0.0299j,  0.4367-0.7640j]]])
        >>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(k,dtype=q.dtype))
        True
    """
    assert k<=n
    assert u.size(-1)==(n**2+n)
    assert (0<=u).all()
    assert (u<=1).all()
    num_tril = n*(n-1)//2
    alpha = agsutil.icdf_std_normal(u[...,n:2*n])*np.sqrt(2)
    beta = torch.complex(
        agsutil.icdf_std_normal(u[...,2*n:2*n+num_tril]),
        agsutil.icdf_std_normal(u[...,2*n+num_tril:]))
    theta = 2*np.pi*u[...,:n]
    il0,il1 = torch.tril_indices(n,n,offset=-1,device=u.device)
    v = torch.zeros((*u.shape[:-1],n,n),dtype=beta.dtype,device=u.device)
    diag_indices = torch.arange(n,device=u.device)
    v[...,diag_indices,diag_indices] = alpha.to(beta.dtype)
    v[...,il0,il1] = beta
    v += v.tril(-1).conj().transpose(dim0=-2,dim1=-1)
    # assert torch.allclose(v[...,torch.arange(n,device=u.device),torch.arange(n,device=u.device)],alpha.to(v.dtype))
    gamma,q = torch.linalg.eigh(v)
    phases = torch.complex(torch.cos(theta),torch.sin(theta))
    q = q*phases[...,None,:]
    # assert torch.allclose(torch.einsum("...ij,...j,...kj->...ik",q,gamma,q.conj()),v)
    return q[...,:k]

tf_Spn_SVD 

tf_Spn_SVD(u, n, k)

Transform uniform samples to random symplectic unitary matrices using SVD.

Parameters:

Name Type Description Default
u Tensor

Uniformly distributed samples in [0, 1] of shape (..., 4*n**2).

 required
n int

Half the number of rows.

 required
k int

Half the number of columns with k<=n.

 required

Returns:

Name Type Description
x Tensor

A batch of random symplectic unitary matrices of shape (..., 2n, 2k).

Examples:

>>> torch.set_default_dtype(torch.float64) 
>>> rng = torch.Generator().manual_seed(7)

>>> delta = torch.arange(1,4,dtype=float)
>>> delta = delta/torch.linalg.norm(delta)
>>> n = delta.size(-1)
>>> q = tf_Spn_SVD(torch.rand((2,4*n**2),generator=rng),n,n)
>>> q.shape
torch.Size([2, 6, 6])
>>> q
tensor([[[-0.3855-0.0066j,  0.2846+0.1340j,  0.0531+0.4453j, -0.0713-0.0289j,
          -0.0862-0.1470j, -0.5170-0.4990j],
         [-0.2564-0.3464j, -0.4935+0.2373j,  0.3129-0.0210j, -0.1699+0.2938j,
          -0.3639+0.4087j, -0.0355-0.0141j],
         [ 0.4523+0.2013j, -0.3856+0.2435j,  0.2998+0.1332j,  0.5145+0.1862j,
           0.0854-0.2397j, -0.1482-0.2306j],
         [ 0.0713-0.0289j,  0.0862-0.1470j,  0.5170-0.4990j, -0.3855+0.0066j,
           0.2846-0.1340j,  0.0531-0.4453j],
         [ 0.1699+0.2938j,  0.3639+0.4087j,  0.0355-0.0141j, -0.2564+0.3464j,
          -0.4935-0.2373j,  0.3129+0.0210j],
         [-0.5145+0.1862j, -0.0854-0.2397j,  0.1482-0.2306j,  0.4523-0.2013j,
          -0.3856-0.2435j,  0.2998-0.1332j]],

        [[-0.4396-0.3859j, -0.1858+0.4638j, -0.1706-0.4717j,  0.1578-0.1175j,
          -0.1740+0.0735j,  0.2780-0.0707j],
         [ 0.0092+0.2218j, -0.1126+0.4596j,  0.3860+0.1825j, -0.4997-0.0507j,
           0.0025+0.3514j,  0.2753+0.3050j],
         [ 0.3743+0.2976j,  0.2388+0.3247j,  0.0849-0.3571j,  0.2850-0.0874j,
          -0.4425+0.0953j, -0.3876+0.1736j],
         [-0.1578-0.1175j,  0.1740+0.0735j, -0.2780-0.0707j, -0.4396+0.3859j,
          -0.1858-0.4638j, -0.1706+0.4717j],
         [ 0.4997-0.0507j, -0.0025+0.3514j, -0.2753+0.3050j,  0.0092-0.2218j,
          -0.1126-0.4596j,  0.3860-0.1825j],
         [-0.2850-0.0874j,  0.4425+0.0953j,  0.3876+0.1736j,  0.3743-0.2976j,
           0.2388-0.3247j,  0.0849+0.3571j]]])
>>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(2*n,dtype=q.dtype))
True

>>> k = 2 
>>> q = tf_Spn_SVD(torch.rand((5,4*n**2),generator=rng),n,k)
>>> q.shape
torch.Size([5, 6, 4])
>>> q
tensor([[[-0.0633-0.1036j, -0.2454-0.0317j, -0.1075+0.2173j,  0.0357-0.2214j],
         [-0.7398+0.2957j, -0.4981-0.1796j, -0.1062+0.0984j,  0.1297+0.0526j],
         [ 0.0691+0.3965j,  0.2245-0.3309j,  0.0134-0.3295j,  0.4265+0.4968j],
         [ 0.1075+0.2173j, -0.0357-0.2214j, -0.0633+0.1036j, -0.2454+0.0317j],
         [ 0.1062+0.0984j, -0.1297+0.0526j, -0.7398-0.2957j, -0.4981+0.1796j],
         [-0.0134-0.3295j, -0.4265+0.4968j,  0.0691-0.3965j,  0.2245+0.3309j]],

        [[-0.0486-0.0043j,  0.1458+0.0475j,  0.2022-0.0570j,  0.0165+0.1342j],
         [ 0.2863-0.1145j, -0.3776+0.5283j, -0.0927-0.4884j,  0.3623-0.2821j],
         [-0.5816+0.0233j, -0.3505+0.1184j, -0.1007+0.5121j,  0.3969-0.1769j],
         [-0.2022-0.0570j, -0.0165+0.1342j, -0.0486+0.0043j,  0.1458-0.0475j],
         [ 0.0927-0.4884j, -0.3623-0.2821j,  0.2863+0.1145j, -0.3776-0.5283j],
         [ 0.1007+0.5121j, -0.3969-0.1769j, -0.5816-0.0233j, -0.3505-0.1184j]],

        [[ 0.1763-0.0398j, -0.5018-0.4715j,  0.1128+0.1180j, -0.0184+0.4873j],
         [-0.6141-0.1638j,  0.0427-0.1202j,  0.1371+0.1876j,  0.1736-0.2845j],
         [ 0.4854-0.2710j, -0.1455-0.1455j,  0.1575-0.3858j, -0.1879-0.2882j],
         [-0.1128+0.1180j,  0.0184+0.4873j,  0.1763+0.0398j, -0.5018+0.4715j],
         [-0.1371+0.1876j, -0.1736-0.2845j, -0.6141+0.1638j,  0.0427+0.1202j],
         [-0.1575-0.3858j,  0.1879-0.2882j,  0.4854+0.2710j, -0.1455+0.1455j]],

        [[ 0.0197-0.1907j, -0.1932-0.7409j, -0.1938-0.1155j,  0.4894-0.0170j],
         [ 0.1002+0.0349j, -0.2719+0.1094j,  0.1073+0.0866j,  0.0376-0.0229j],
         [-0.2872+0.4903j, -0.0111+0.1788j, -0.7311-0.1571j,  0.1123+0.2034j],
         [ 0.1938-0.1155j, -0.4894-0.0170j,  0.0197+0.1907j, -0.1932+0.7409j],
         [-0.1073+0.0866j, -0.0376-0.0229j,  0.1002-0.0349j, -0.2719-0.1094j],
         [ 0.7311-0.1571j, -0.1123+0.2034j, -0.2872-0.4903j, -0.0111-0.1788j]],

        [[ 0.1997+0.0214j, -0.4774+0.0789j,  0.3295-0.2111j,  0.1172-0.4257j],
         [ 0.0496+0.7202j, -0.1458-0.1148j, -0.1592+0.0259j, -0.4487-0.0895j],
         [-0.1664+0.2085j, -0.1348+0.0409j, -0.0706-0.4281j,  0.5542+0.0138j],
         [-0.3295-0.2111j, -0.1172-0.4257j,  0.1997-0.0214j, -0.4774-0.0789j],
         [ 0.1592+0.0259j,  0.4487-0.0895j,  0.0496-0.7202j, -0.1458+0.1148j],
         [ 0.0706-0.4281j, -0.5542+0.0138j, -0.1664-0.2085j, -0.1348-0.0409j]]])
>>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(2*k,dtype=q.dtype))
True
Source code in rmatqmcpy/tf_On_Un_Spn.py 
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
def tf_Spn_SVD(u, n, k):
    r"""
    Transform uniform samples to random symplectic unitary matrices using SVD.

    Args:
        u (torch.Tensor): Uniformly distributed samples in [0, 1] of shape `(..., 4*n**2)`.
        n (int): Half the number of rows.
        k (int): Half the number of columns with `k<=n`.

    Returns:
        x (torch.Tensor): A batch of random symplectic unitary matrices of shape `(..., 2n, 2k)`.

    Examples:
        >>> torch.set_default_dtype(torch.float64) 
        >>> rng = torch.Generator().manual_seed(7)

        >>> delta = torch.arange(1,4,dtype=float)
        >>> delta = delta/torch.linalg.norm(delta)
        >>> n = delta.size(-1)
        >>> q = tf_Spn_SVD(torch.rand((2,4*n**2),generator=rng),n,n)
        >>> q.shape
        torch.Size([2, 6, 6])
        >>> q
        tensor([[[-0.3855-0.0066j,  0.2846+0.1340j,  0.0531+0.4453j, -0.0713-0.0289j,
                  -0.0862-0.1470j, -0.5170-0.4990j],
                 [-0.2564-0.3464j, -0.4935+0.2373j,  0.3129-0.0210j, -0.1699+0.2938j,
                  -0.3639+0.4087j, -0.0355-0.0141j],
                 [ 0.4523+0.2013j, -0.3856+0.2435j,  0.2998+0.1332j,  0.5145+0.1862j,
                   0.0854-0.2397j, -0.1482-0.2306j],
                 [ 0.0713-0.0289j,  0.0862-0.1470j,  0.5170-0.4990j, -0.3855+0.0066j,
                   0.2846-0.1340j,  0.0531-0.4453j],
                 [ 0.1699+0.2938j,  0.3639+0.4087j,  0.0355-0.0141j, -0.2564+0.3464j,
                  -0.4935-0.2373j,  0.3129+0.0210j],
                 [-0.5145+0.1862j, -0.0854-0.2397j,  0.1482-0.2306j,  0.4523-0.2013j,
                  -0.3856-0.2435j,  0.2998-0.1332j]],
        <BLANKLINE>
                [[-0.4396-0.3859j, -0.1858+0.4638j, -0.1706-0.4717j,  0.1578-0.1175j,
                  -0.1740+0.0735j,  0.2780-0.0707j],
                 [ 0.0092+0.2218j, -0.1126+0.4596j,  0.3860+0.1825j, -0.4997-0.0507j,
                   0.0025+0.3514j,  0.2753+0.3050j],
                 [ 0.3743+0.2976j,  0.2388+0.3247j,  0.0849-0.3571j,  0.2850-0.0874j,
                  -0.4425+0.0953j, -0.3876+0.1736j],
                 [-0.1578-0.1175j,  0.1740+0.0735j, -0.2780-0.0707j, -0.4396+0.3859j,
                  -0.1858-0.4638j, -0.1706+0.4717j],
                 [ 0.4997-0.0507j, -0.0025+0.3514j, -0.2753+0.3050j,  0.0092-0.2218j,
                  -0.1126-0.4596j,  0.3860-0.1825j],
                 [-0.2850-0.0874j,  0.4425+0.0953j,  0.3876+0.1736j,  0.3743-0.2976j,
                   0.2388-0.3247j,  0.0849+0.3571j]]])
        >>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(2*n,dtype=q.dtype))
        True

        >>> k = 2 
        >>> q = tf_Spn_SVD(torch.rand((5,4*n**2),generator=rng),n,k)
        >>> q.shape
        torch.Size([5, 6, 4])
        >>> q
        tensor([[[-0.0633-0.1036j, -0.2454-0.0317j, -0.1075+0.2173j,  0.0357-0.2214j],
                 [-0.7398+0.2957j, -0.4981-0.1796j, -0.1062+0.0984j,  0.1297+0.0526j],
                 [ 0.0691+0.3965j,  0.2245-0.3309j,  0.0134-0.3295j,  0.4265+0.4968j],
                 [ 0.1075+0.2173j, -0.0357-0.2214j, -0.0633+0.1036j, -0.2454+0.0317j],
                 [ 0.1062+0.0984j, -0.1297+0.0526j, -0.7398-0.2957j, -0.4981+0.1796j],
                 [-0.0134-0.3295j, -0.4265+0.4968j,  0.0691-0.3965j,  0.2245+0.3309j]],
        <BLANKLINE>
                [[-0.0486-0.0043j,  0.1458+0.0475j,  0.2022-0.0570j,  0.0165+0.1342j],
                 [ 0.2863-0.1145j, -0.3776+0.5283j, -0.0927-0.4884j,  0.3623-0.2821j],
                 [-0.5816+0.0233j, -0.3505+0.1184j, -0.1007+0.5121j,  0.3969-0.1769j],
                 [-0.2022-0.0570j, -0.0165+0.1342j, -0.0486+0.0043j,  0.1458-0.0475j],
                 [ 0.0927-0.4884j, -0.3623-0.2821j,  0.2863+0.1145j, -0.3776-0.5283j],
                 [ 0.1007+0.5121j, -0.3969-0.1769j, -0.5816-0.0233j, -0.3505-0.1184j]],
        <BLANKLINE>
                [[ 0.1763-0.0398j, -0.5018-0.4715j,  0.1128+0.1180j, -0.0184+0.4873j],
                 [-0.6141-0.1638j,  0.0427-0.1202j,  0.1371+0.1876j,  0.1736-0.2845j],
                 [ 0.4854-0.2710j, -0.1455-0.1455j,  0.1575-0.3858j, -0.1879-0.2882j],
                 [-0.1128+0.1180j,  0.0184+0.4873j,  0.1763+0.0398j, -0.5018+0.4715j],
                 [-0.1371+0.1876j, -0.1736-0.2845j, -0.6141+0.1638j,  0.0427+0.1202j],
                 [-0.1575-0.3858j,  0.1879-0.2882j,  0.4854+0.2710j, -0.1455+0.1455j]],
        <BLANKLINE>
                [[ 0.0197-0.1907j, -0.1932-0.7409j, -0.1938-0.1155j,  0.4894-0.0170j],
                 [ 0.1002+0.0349j, -0.2719+0.1094j,  0.1073+0.0866j,  0.0376-0.0229j],
                 [-0.2872+0.4903j, -0.0111+0.1788j, -0.7311-0.1571j,  0.1123+0.2034j],
                 [ 0.1938-0.1155j, -0.4894-0.0170j,  0.0197+0.1907j, -0.1932+0.7409j],
                 [-0.1073+0.0866j, -0.0376-0.0229j,  0.1002-0.0349j, -0.2719-0.1094j],
                 [ 0.7311-0.1571j, -0.1123+0.2034j, -0.2872-0.4903j, -0.0111-0.1788j]],
        <BLANKLINE>
                [[ 0.1997+0.0214j, -0.4774+0.0789j,  0.3295-0.2111j,  0.1172-0.4257j],
                 [ 0.0496+0.7202j, -0.1458-0.1148j, -0.1592+0.0259j, -0.4487-0.0895j],
                 [-0.1664+0.2085j, -0.1348+0.0409j, -0.0706-0.4281j,  0.5542+0.0138j],
                 [-0.3295-0.2111j, -0.1172-0.4257j,  0.1997-0.0214j, -0.4774-0.0789j],
                 [ 0.1592+0.0259j,  0.4487-0.0895j,  0.0496-0.7202j, -0.1458+0.1148j],
                 [ 0.0706-0.4281j, -0.5542+0.0138j, -0.1664-0.2085j, -0.1348-0.0409j]]])
        >>> torch.allclose(torch.einsum("...ji,...jk->...ik",q,q.conj()),torch.eye(2*k,dtype=q.dtype))
        True

    """ 
    assert k<=n
    assert u.size(-1)==(4*n**2)
    assert (0<=u).all()
    assert (u<=1).all()
    u = u.reshape((*u.shape[:-1],n,n,4))
    v = agsutil.icdf_std_normal(u)
    B = torch.complex(v[..., 0],v[..., 1])
    C = torch.complex(v[..., 2],v[..., 3])
    v_symplectic = torch.cat([
        torch.cat([B,              C            ], dim=-1),
        torch.cat([-torch.conj(C), torch.conj(B)], dim=-1)], dim=-2)
    U,_,Vh = torch.linalg.svd(v_symplectic)
    q = torch.einsum("...ij,...jk->...ik",U,Vh)
    ivec = torch.cat([torch.arange(k),torch.arange(n,n+k)],dim=0)
    return q[...,ivec]