Skip to content

perm

The quotonic.perm module includes functions and classes used to compute matrix permanents as efficiently as possible when using jax as a backend. Currently, there is support both for Ryser's algorithm and the Balasubramanian-Bax-Franklin-Glynn (BBFF) algorithm, each using Gray code ordering.

This code was inspired by Piquasso, The Walrus, and Cascaded Optical Systems Approach to Neural Networks (CasOptAx) as originally designed for use in J. R. Basani et al., "Universal logical quantum photonic neural network processor via cavity-assisted interactions", npj Quantum Inf 11, 142 (2025).

prep_gray_code(i)

Preparation of Gray code for computing matrix permanents.

This function is wrapped with jax.vmap such that it can be used vectorially.

Parameters:

Name Type Description Default
i int

index typically used in loops that prepare Gray code

 required

Returns:

Name Type Description
gray_diff int

the difference between the old Gray value and the new
direction int

the direction of the algorithm, either +1 or -1
Source code in src/quotonic/perm.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
@vmap
def prep_gray_code(i: int) -> tuple[int, int]:
    """Preparation of Gray code for computing matrix permanents.

    This function is wrapped with `jax.vmap` such that it can be used vectorially.

    Args:
        i: index typically used in loops that prepare Gray code

    Returns:
        gray_diff: the difference between the old Gray value and the new
        direction: the direction of the algorithm, either +1 or -1
    """
    old_gray = i ^ (i // 2)
    new_gray = (i + 1) ^ ((i + 1) // 2)
    gray_diff = old_gray ^ new_gray
    direction = lax.cond(
        old_gray > new_gray,
        lambda: 1,
        lambda: -1,
    )
    return gray_diff, direction

calc_perm_ryser(U)

Compute the permanent of a square matrix \(\mathbf{U}\) using Ryser's algorithm with Gray code ordering.

This function is wrapped with jax.jit such that it can be compiled at runtime.

Parameters:

Name Type Description Default
U jnp_ndarray

square matrix whose permanent is to be computed

 required

Returns:

Name Type Description
perm DTypeLike

permanent of the given square matrix
Source code in src/quotonic/perm.py 
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
@jit
def calc_perm_ryser(U: jnp_ndarray) -> DTypeLike:
    """Compute the permanent of a square matrix $\\mathbf{U}$ using Ryser's algorithm with Gray code ordering.

    This function is wrapped with `jax.jit` such that it can be compiled at runtime.

    Args:
        U: square matrix whose permanent is to be computed

    Returns:
        perm: permanent of the given square matrix
    """

    n = U.shape[0]
    two_to_n = 2**n

    # prepare Gray code for the permanent calculation
    gray_diff, direction = prep_gray_code(jnp.arange(two_to_n - 1, dtype=jnp.int16))  # type: ignore
    gray_diff_ind = jnp.array(jnp.log2(gray_diff), dtype=jnp.int16)
    sign = jnp.resize(jnp.array([-1, 1], dtype=jnp.int16), two_to_n - 1)

    # calculate permanent by vectorizing Ryser's algorithm
    perm: DTypeLike = jnp.sum(
        sign * jnp.prod(jnp.cumsum(vmap(lambda ind, direc: U[ind] * direc)(gray_diff_ind, direction), axis=0), axis=1)
    )
    return perm

calc_perm_bbfg(U)

Compute the permanent of a square matrix \(\mathbf{U}\) using the BBFG algorithm with Gray code ordering.

This function is wrapped with jax.jit such that it can be compiled at runtime.

Parameters:

Name Type Description Default
U jnp_ndarray

square matrix whose permanent is to be computed

 required

Returns:

Name Type Description
perm DTypeLike

permanent of the given square matrix
Source code in src/quotonic/perm.py 
 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
@jit
def calc_perm_bbfg(U: jnp_ndarray) -> DTypeLike:
    """Compute the permanent of a square matrix $\\mathbf{U}$ using the BBFG algorithm with Gray code ordering.

    This function is wrapped with `jax.jit` such that it can be compiled at runtime.

    Args:
        U: square matrix whose permanent is to be computed

    Returns:
        perm: permanent of the given square matrix
    """

    n = U.shape[0]
    N = 2 ** (n - 1)

    # prepare Gray code for the permanent calculation
    gray_diff, direction = prep_gray_code(jnp.arange(N - 1, dtype=jnp.int16))  # type: ignore
    gray_diff_ind = jnp.array(jnp.log2(gray_diff), dtype=jnp.int16)
    sign = jnp.resize(jnp.array([-1, 1], dtype=jnp.int16), N - 1)

    # calculate permanent by vectorizing the BBFG algorithm
    perm: DTypeLike = (
        jnp.prod(jnp.sum(U, axis=0))
        + jnp.sum(
            sign
            * (
                jnp.prod(
                    jnp.sum(U, axis=0)
                    + jnp.cumsum(
                        vmap(lambda ind, direc: U[ind] * direc * 2)(gray_diff_ind, direction),
                        axis=0,
                    ),
                    axis=1,
                )
            )
        )
    ) / N
    return perm

calc_perm(U, algo='bbfg')

Compute the permanent of a square matrix \(\mathbf{U}\).

Parameters:

Name Type Description Default
U jnp_ndarray

square matrix whose permanent is to be computed

 required
algo str

algorithm to compute the permanent with if the matrix dim is greater than 3, either "bbfg" or "ryser"

 'bbfg'

Returns:

Name Type Description
perm DTypeLike

permanent of the given square matrix
Source code in src/quotonic/perm.py 
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
@partial(jit, static_argnums=(1,))
def calc_perm(U: jnp_ndarray, algo: str = "bbfg") -> DTypeLike:
    """Compute the permanent of a square matrix $\\mathbf{U}$.

    Args:
        U: square matrix whose permanent is to be computed
        algo: algorithm to compute the permanent with if the matrix dim is greater than 3, either "bbfg" or "ryser"

    Returns:
        perm: permanent of the given square matrix
    """

    # extract the dimension of the square matrix U
    Ushape = jnp.shape(U)
    assert Ushape[0] == Ushape[1], "Matrix must be square"
    assert Ushape[0] > 0, "Matrix must have elements"
    N = Ushape[0]

    if N == 1:
        return U[0, 0]

    if N == 2:
        return U[0, 0] * U[1, 1] + U[0, 1] * U[1, 0]

    if N == 3:
        return (
            U[0, 2] * U[1, 1] * U[2, 0]
            + U[0, 1] * U[1, 2] * U[2, 0]
            + U[0, 2] * U[1, 0] * U[2, 1]
            + U[0, 0] * U[1, 2] * U[2, 1]
            + U[0, 1] * U[1, 0] * U[2, 2]
            + U[0, 0] * U[1, 1] * U[2, 2]
        )

    perm: DTypeLike = calc_perm_bbfg(U) if algo == "bbfg" else calc_perm_ryser(U)
    return perm

Permanent

Wrapper class for computing permanents of matrices of constant dimension \(n\) while the Gray code overhead is stored in memory.

Attributes:

Name Type Description
n int

dimension of the square matrices, \(n\)
perm callable

function that computes the permanent of a given \(n\times n\) matrix \(\mathbf{U}\)
gray_diff_ind jnp_ndarray

array of matrix indices for computing permanents using Gray code ordering, defaults to an empty array if \(n < 3\)
direction jnp_ndarray

array of factors to apply in individual steps of the permanent calculation algorithms, defaults to an empty array if \(n < 3\)
sign jnp_ndarray

array of \(\pm 1\) factors to apply to the results of individual steps of the permanent calculation algorithms, defaults to an empty array if \(n < 3\)
N int

if the BBFG algorithm is selected, \(N = 2^{n-1}\), otherwise, it is unused and defaults to 0
Source code in src/quotonic/perm.py 
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
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
class Permanent:
    """Wrapper class for computing permanents of matrices of constant dimension $n$ while the Gray code overhead is
    stored in memory.

    Attributes:
        n (int): dimension of the square matrices, $n$
        perm (callable): function that computes the permanent of a given $n\\times n$ matrix $\\mathbf{U}$
        gray_diff_ind (jnp_ndarray): array of matrix indices for computing permanents using Gray code ordering,
            defaults to an empty array if $n < 3$
        direction (jnp_ndarray): array of factors to apply in individual steps of the permanent calculation
            algorithms, defaults to an empty array if $n < 3$
        sign (jnp_ndarray): array of $\\pm 1$ factors to apply to the results of individual steps of the permanent
            calculation algorithms, defaults to an empty array if $n < 3$
        N (int): if the BBFG algorithm is selected, $N = 2^{n-1}$, otherwise, it is unused and defaults to 0
    """

    def __init__(self, n: int, algo: str = "bbfg") -> None:
        """Initialization of a Permanent instance.

        Args:
            n: dimension of the square matrices whose permanents are to be computed, $n$
            algo: algorithm to compute the permanent with if the matrix dimension is greater than 3, either "bbfg" or "ryser"
        """

        # check the validity of the provided arguments
        assert n > 0, "Matrices must have elements to compute a permanent"
        assert (algo == "bbfg") or (algo == "ryser"), "The only algorithm options are 'bbfg' or 'ryser'"

        # store the dimension of the square matrices, then compute overhead if necessary
        self.n = n
        if n < 3:
            self.gray_diff_ind = jnp.array(())
            self.direction = jnp.array(())
            self.sign = jnp.array(())
            self.N = 0
            if n == 1:
                self.perm = jit(lambda U: U[0, 0])
            elif n == 2:
                self.perm = jit(lambda U: U[0, 0] * U[1, 1] + U[0, 1] * U[1, 0])
            elif n == 3:
                self.perm = jit(
                    lambda U: U[0, 2] * U[1, 1] * U[2, 0]
                    + U[0, 1] * U[1, 2] * U[2, 0]
                    + U[0, 2] * U[1, 0] * U[2, 1]
                    + U[0, 0] * U[1, 2] * U[2, 1]
                    + U[0, 1] * U[1, 0] * U[2, 2]
                    + U[0, 0] * U[1, 1] * U[2, 2]
                )
        else:
            if algo == "bbfg":
                self.prep_gray_code_bbfg()
                self.perm = self.perm_bbfg
            elif algo == "ryser":
                self.prep_gray_code_ryser()
                self.perm = self.perm_ryser
                self.N = 0

    def prep_gray_code_bbfg(self) -> None:
        """Preparation of Gray code for computing matrix permanents using the BBFG algorithm."""
        self.N = 2 ** (self.n - 1)
        gray_diff, direction = prep_gray_code(jnp.arange(self.N - 1, dtype=jnp.int16))  # type: ignore
        self.direction = 2 * direction  # type: ignore
        self.gray_diff_ind = jnp.array(jnp.log2(gray_diff), dtype=jnp.int16)
        self.sign = jnp.resize(jnp.array([-1, 1], dtype=jnp.int16), (self.N - 1,))

    def prep_gray_code_ryser(self) -> None:
        """Preparation of Gray code for computing matrix permanents using Ryser's algorithm."""
        two_to_n = 2**self.n
        gray_diff, self.direction = prep_gray_code(jnp.arange(two_to_n - 1, dtype=jnp.int16))  # type: ignore
        self.gray_diff_ind = jnp.array(jnp.log2(gray_diff), dtype=jnp.int16)
        self.sign = jnp.resize(jnp.array([-1, 1], dtype=jnp.int16), (two_to_n - 1,))

    @partial(jit, static_argnums=(0,))
    def perm_bbfg(self, U: jnp_ndarray) -> DTypeLike:
        """Compute the permanent of a square matrix $\\mathbf{U}$ using the BBFG algorithm with Gray code ordering.

        Args:
            U: square matrix whose permanent is to be computed

        Returns:
            perm: permanent of the given square matrix
        """
        return (
            jnp.prod(jnp.sum(U, axis=0))
            + jnp.sum(
                self.sign
                * (
                    jnp.prod(
                        jnp.sum(U, axis=0)
                        + jnp.cumsum(
                            vmap(lambda ind, direc: U[ind] * direc)(self.gray_diff_ind, self.direction),
                            axis=0,
                        ),
                        axis=1,
                    )
                )
            )
        ) / self.N

    @partial(jit, static_argnums=(0,))
    def perm_ryser(self, U: jnp_ndarray) -> DTypeLike:
        """Compute the permanent of a square matrix $\\mathbf{U}$ using Ryser's algorithm with Gray code ordering.

        Args:
            U: square matrix whose permanent is to be computed

        Returns:
            perm: permanent of the given square matrix
        """
        return jnp.sum(
            self.sign
            * jnp.prod(
                jnp.cumsum(vmap(lambda ind, direc: U[ind] * direc)(self.gray_diff_ind, self.direction), axis=0), axis=1
            )
        )

__init__(n, algo='bbfg')

Initialization of a Permanent instance.

Parameters:

Name Type Description Default
n int

dimension of the square matrices whose permanents are to be computed, \(n\)

 required
algo str

algorithm to compute the permanent with if the matrix dimension is greater than 3, either "bbfg" or "ryser"

 'bbfg'
Source code in src/quotonic/perm.py 
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
def __init__(self, n: int, algo: str = "bbfg") -> None:
    """Initialization of a Permanent instance.

    Args:
        n: dimension of the square matrices whose permanents are to be computed, $n$
        algo: algorithm to compute the permanent with if the matrix dimension is greater than 3, either "bbfg" or "ryser"
    """

    # check the validity of the provided arguments
    assert n > 0, "Matrices must have elements to compute a permanent"
    assert (algo == "bbfg") or (algo == "ryser"), "The only algorithm options are 'bbfg' or 'ryser'"

    # store the dimension of the square matrices, then compute overhead if necessary
    self.n = n
    if n < 3:
        self.gray_diff_ind = jnp.array(())
        self.direction = jnp.array(())
        self.sign = jnp.array(())
        self.N = 0
        if n == 1:
            self.perm = jit(lambda U: U[0, 0])
        elif n == 2:
            self.perm = jit(lambda U: U[0, 0] * U[1, 1] + U[0, 1] * U[1, 0])
        elif n == 3:
            self.perm = jit(
                lambda U: U[0, 2] * U[1, 1] * U[2, 0]
                + U[0, 1] * U[1, 2] * U[2, 0]
                + U[0, 2] * U[1, 0] * U[2, 1]
                + U[0, 0] * U[1, 2] * U[2, 1]
                + U[0, 1] * U[1, 0] * U[2, 2]
                + U[0, 0] * U[1, 1] * U[2, 2]
            )
    else:
        if algo == "bbfg":
            self.prep_gray_code_bbfg()
            self.perm = self.perm_bbfg
        elif algo == "ryser":
            self.prep_gray_code_ryser()
            self.perm = self.perm_ryser
            self.N = 0

prep_gray_code_bbfg()

Preparation of Gray code for computing matrix permanents using the BBFG algorithm.

Source code in src/quotonic/perm.py 
210
211
212
213
214
215
216
def prep_gray_code_bbfg(self) -> None:
    """Preparation of Gray code for computing matrix permanents using the BBFG algorithm."""
    self.N = 2 ** (self.n - 1)
    gray_diff, direction = prep_gray_code(jnp.arange(self.N - 1, dtype=jnp.int16))  # type: ignore
    self.direction = 2 * direction  # type: ignore
    self.gray_diff_ind = jnp.array(jnp.log2(gray_diff), dtype=jnp.int16)
    self.sign = jnp.resize(jnp.array([-1, 1], dtype=jnp.int16), (self.N - 1,))

prep_gray_code_ryser()

Preparation of Gray code for computing matrix permanents using Ryser's algorithm.

Source code in src/quotonic/perm.py 
218
219
220
221
222
223
def prep_gray_code_ryser(self) -> None:
    """Preparation of Gray code for computing matrix permanents using Ryser's algorithm."""
    two_to_n = 2**self.n
    gray_diff, self.direction = prep_gray_code(jnp.arange(two_to_n - 1, dtype=jnp.int16))  # type: ignore
    self.gray_diff_ind = jnp.array(jnp.log2(gray_diff), dtype=jnp.int16)
    self.sign = jnp.resize(jnp.array([-1, 1], dtype=jnp.int16), (two_to_n - 1,))

perm_bbfg(U)

Compute the permanent of a square matrix \(\mathbf{U}\) using the BBFG algorithm with Gray code ordering.

Parameters:

Name Type Description Default
U jnp_ndarray

square matrix whose permanent is to be computed

 required

Returns:

Name Type Description
perm DTypeLike

permanent of the given square matrix
Source code in src/quotonic/perm.py 
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
@partial(jit, static_argnums=(0,))
def perm_bbfg(self, U: jnp_ndarray) -> DTypeLike:
    """Compute the permanent of a square matrix $\\mathbf{U}$ using the BBFG algorithm with Gray code ordering.

    Args:
        U: square matrix whose permanent is to be computed

    Returns:
        perm: permanent of the given square matrix
    """
    return (
        jnp.prod(jnp.sum(U, axis=0))
        + jnp.sum(
            self.sign
            * (
                jnp.prod(
                    jnp.sum(U, axis=0)
                    + jnp.cumsum(
                        vmap(lambda ind, direc: U[ind] * direc)(self.gray_diff_ind, self.direction),
                        axis=0,
                    ),
                    axis=1,
                )
            )
        )
    ) / self.N

perm_ryser(U)

Compute the permanent of a square matrix \(\mathbf{U}\) using Ryser's algorithm with Gray code ordering.

Parameters:

Name Type Description Default
U jnp_ndarray

square matrix whose permanent is to be computed

 required

Returns:

Name Type Description
perm DTypeLike

permanent of the given square matrix
Source code in src/quotonic/perm.py 
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
@partial(jit, static_argnums=(0,))
def perm_ryser(self, U: jnp_ndarray) -> DTypeLike:
    """Compute the permanent of a square matrix $\\mathbf{U}$ using Ryser's algorithm with Gray code ordering.

    Args:
        U: square matrix whose permanent is to be computed

    Returns:
        perm: permanent of the given square matrix
    """
    return jnp.sum(
        self.sign
        * jnp.prod(
            jnp.cumsum(vmap(lambda ind, direc: U[ind] * direc)(self.gray_diff_ind, self.direction), axis=0), axis=1
        )
    )

EmptyPermanent

Placeholder permanent class for mypy typing.

Source code in src/quotonic/perm.py 
270
271
272
273
274
275
276
277
278
279
class EmptyPermanent:
    """Placeholder permanent class for mypy typing."""

    def __init__(self) -> None:
        """Initialization of an EmptyPermanent instance."""
        self.fake = True

    def perm(self, U: jnp_ndarray) -> jnp_ndarray:
        """Fake method required for the placeholder."""
        return U[0, 0] if self.fake else U[0, 0]

__init__()

Initialization of an EmptyPermanent instance.

Source code in src/quotonic/perm.py 
273
274
275
def __init__(self) -> None:
    """Initialization of an EmptyPermanent instance."""
    self.fake = True

perm(U)

Fake method required for the placeholder.

Source code in src/quotonic/perm.py 
277
278
279
def perm(self, U: jnp_ndarray) -> jnp_ndarray:
    """Fake method required for the placeholder."""
    return U[0, 0] if self.fake else U[0, 0]