Technical Field

[0001] The present disclosure relates to methods and systems for addressing binary optimization problems. More particularly, the present disclosure relates to methods and systems that utilize one or more boson sampling devices to determine solutions to binary optimization problems.

Background

[0002] Binary optimization problems are a subclass of combinatorial optimization problems in which the variables are restricted to one of two values. A binary optimization problem can generally be stated as a problem of reducing (e.g., minimizing) or increasing (e.g., maximizing) an objective function of variables that can take one of two values, for example 0 or 1.

[0003] One famous binary optimization problem is the Quadratic Unconstrained Binary Optimization (QUBO) problem, which is also known as the Unconstrained Binary Quadratic Programming (UBQP) problem. The task is to find a binary sequence b of length L that minimises the objective function b ^T Qb where Q is an L X L symmetric matrix having real values. Many problems of real-world relevance can be cast into a QUBO form. QUBO is known to be an NP hard problem and so is intractable on a classical computer when the number L of binary variables is large.

[0004] Systems and methods for finding solutions to binary optimization problems are desirable.

Summary

[0005] Some embodiments relate to a system for determining a solution to a binary optimization problem. As used herein, “solution” may refer to a value that exactly minimizes (or maximizes) an objective function (e.g., a globally optimal solution). However, the term should not be limited to this definition. For example, a “solution” may refer to a value that reduces (or increases) an objective function relative to an initial value. In another example, a “solution” results in a local minimum or maximum of the objective function.

[0006] A binary sequence is understood to mean any ordered sequence of elements that can take one of two values. In some examples, elements of a binary sequence may take a value of 0 or 1. In other examples, elements of a binary sequence may take a value of —1 or 1. [0007] The system comprises a boson sampler and one or more classical processors.

[0008] The boson sampler comprises a set of one or more single photon sources, the set configured to generate an input photonic state comprising a plurality of input modes. The boson sampler further comprises a configurable interferometer for generating an output photonic state comprising a plurality of output modes by interfering photons of the input photonic state. The boson sampler further comprises an arrangement of one or more photon detectors, the arrangement configured to produce a measurement outcome indicative of a presence or absence of photons in output modes of the output photonic state.

[0009] The one or more processors are configured to (i) initialize a set of parameter values, the set of parameter values for defining the configuration of the interferometer. The one or more processors are further configured to (ii) for at least one selected parameter, use the boson sampler to determine a gradient of a cost function with respect to that selected parameter. The one or more processors are further configured to, using the at least one determined gradient, update the set of parameter values. The one or more processors are further configured to (iv) repeat (ii) and (iii) until a stopping condition is satisfied. The one or more processors are further configured to (v) after the stopping condition is satisfied, control the boson sampler to produce a plurality of measurement outcomes. The one or more processors are further configured to (vi) map each distinct measurement outcome to a binary sequence of a plurality of binary sequences, wherein each element of a binary sequence corresponds to a measured output mode of the output photonic state. The one or more processors are further configured to (vii) evaluate an objective function using binary sequences to which distinct measurement outcomes are mapped to identify a solution to the binary optimization problem.

[0010] Advantageously, the system described herein exploits both a classical computing resource and a quantum computing resource in the form of a boson sampler, to determine a solution to a binary optimization problem. By exploiting bosonic statistics and mapping measurement outcomes to binary sequences, the system is able to handle even large binary optimization problems (those for which the candidate solutions have a large number of elements). Furthermore, by mapping measurement outcomes to binary sequences based on a parity of the number of photons in an output mode or based solely on the presence or absence of photons in an output mode, the boson sampler may be a shallow (i.e. not full-depth) boson sampler. If the input photonic state and the output photonic state each comprise M modes, then the set of parameter values for defining the parametrised unitary transformation of the boson sampler may comprise as few as M — 1 parameter values. As a shallow boson sampler may be utilised, the number of components and therefore the amount of error in the computation can be greatly reduced.

[0011] The one or more processors may comprise one or more central processing units (CPUs), and/or one or more graphics processing units (GPUs) and/or application specific circuitry such as a configured field programmable gate array (FPGA) or application specific integrated circuit (ASIC).

[0012] According to some examples, the arrangement of one or more photon detectors comprises an arrangement of one or more photon number resolving (PNR) detectors. The use of PNR detectors advantageously enables a greater selection of mappings to be used to suit the binary optimization problem being addressed.

[0013] For example, when PNR detectors are used, each element of a binary sequence to which a measurement outcome is mapped may have a value based on a parity of the number of photons in the measured output mode to which that element corresponds .

[0014] When parity mappings are to be used, mapping each measurement outcome to a binary sequence may comprise mapping each measurement outcome to a binary sequence according to a first mapping under which each element of a binary sequence has a first value if the number of photons in the corresponding output mode is odd and a second binary value if the number of photons in the corresponding output mode is even. The one or more processors may be further configured to repeat (i) to (vii) with the complementary parity mapping. In other words, the one or more processors may be further configured to repeat (i) to (vii), wherein mapping each measurement outcome to a binary sequence comprises mapping each measurement outcome to a binary sequence according to a second mapping under which each element of a binary sequence has the first binary value if the number of photons in the corresponding output mode is even and the second binary value if the number of photons in the corresponding output mode is odd. Advantageously, the globally optimal solution to the binary optimization problem is then guaranteed to be in the solution set accessible by the system.

[0015] The set of one or more single-photon sources may be configured to generate a first input photonic state comprising M input modes, each input mode comprising a photon, and the set of one or more single-photon sources may be further configured to generate a second input photonic state comprising M input modes, with M - 1 input modes comprising a photon. The one or more processors may be configured to perform (ii) and (v) with the set of one or more single-photon sources configured to generate one of the first or second input photonic state; and configured to repeat (i) to (vii), wherein (ii) and (v) are performed with the set of one or more single-photon sources configured to generate the other of the first or second input photonic state. In some examples, by using different input photonic states, better convergence on a solution can advantageously be achieved. For example, as will be described further below herein, when parity mappings are used some binary sequences are mapped to from a greater number of measurement outcomes than other binary sequences, and so it is advantageous to perform the method for two different parity mappings and two different input states in order to mitigate such imbalances in the parity mapping.

[0016] In some examples, each element of a binary sequence to which a measurement outcome is mapped may have a value based on whether one or more photons were present or absent in the measured output mode to which that element corresponds. This mapping may be used if the boson sampler comprises PNR detectors. However, an advantage of such a mapping is that number resolving detectors are not required: the boson sampler does not require photon number resolving detectors as only the presence or absence of bosons need be determined, and this means that less sensitive detectors are required. In some examples, the arrangement of one or more photon detectors may comprise an arrangement of one or more threshold detectors, also known as on-off detectors. For example, a threshold detector may comprise an avalanche photodiode detector. When presence/absence mappings are to be used, mapping each measurement outcome to a binary sequence may comprise mapping each measurement outcome to a binary sequence according to a first mapping under which each element of a binary sequence has a first value if one or more photons were present in the corresponding measured output mode and a second value if no photons were present in the corresponding measured output mode. In some examples, the one or more processors are further configured to repeat (i) to (vii) with the complementary absence/presence mapping. In other words, the one or more processors may be further configured to repeat (i) to (vii) wherein mapping each measurement outcome to a binary sequence may comprise mapping each measurement outcome to a binary sequence according to a second mapping under which each element of a binary sequence has the second value if one or more photons were present in the corresponding measured output mode and the first value if no photons were present in the corresponding measured output mode.

[0017] The interferometer may be a shallow interferometer. A shallow interferometer is understood to mean an interferometer that is not a full-depth interferometer. A full-depth interferometer is an interferometer that comprises suitable parametrised multimodal operations such that a photon in any input mode of the input photonic state may be scattered to any output mode of the output photonic state. In a shallow circuit, there are limitations on how the photons of the input modes can be scattered to the output modes of the output photonic state. [0018] The use of a mapping based on the presence or absence of photons or the parity of the number of photons enables a shallow interferometer to be used in the boson sampler. A shallow bosonic circuit uses (e.g., requires) fewer multimodal operations than a full -depth circuit and accordingly the (shallow) boson sampler requires fewer components (e.g. reconfigurable beam splitters) to perform the tasks described herein than a full-depth boson sampler. This provides several advantages: firstly, loss is a major source of error in photonic operations (a key platform for boson sampling) and the fewer components required, the lower the loss; and second, the fewer the components, the smaller the parameter space that is traversed in the process described herein, which means that there are fewer parameter values to update.

[0019] The input photonic state may comprise M input modes (where M is an integer). The interferometer may comprise fewer than — 1) /2 parametrised multimodal operations.

[0020] In some examples, the interferometer may comprise a time bin interferometer. Accordingly, the plurality of input modes may comprise a plurality of temporal modes. For example, the input photonic states may comprise a train of single photons, each spaced apart in time by a fixed amount, such that each photon can be said to be in its own time bin / temporal mode.

[0021] A time bin interferometer may comprise one or more temporal mode coupling devices. In some examples, a temporal mode coupling device may comprise a reconfigurable beam splitter and a delay line, the delay line configured to connect one input port of the reconfigurable beam splitter with one output port of the reconfigurable beam splitter. In some examples, a temporal mode coupling device may comprise a quantum memory. The quantum memory may be configured to place the photons of a temporal mode in a quantum superposition of being stored in the memory and not stored in the memory. The quantum memory may then interact the stored part with a subsequent temporal mode. That is, the quantum memory may be treated provide similar functionality to a reconfigurable beam splitter and delay line.

[0022] In some examples, the interferometer may comprise a spatial interferometer. Accordingly, the plurality of input modes may comprise a plurality of spatial modes. For example, input ports/paths to an interferometer may be thought of as spatial modes.

[0023] A spatial interferometer may comprise M input ports for inputting the M input modes of the input multimodal bosonic state into the interferometer. The interferometer may further comprise M output ports for outputting the M output modes of the output multimodal bosonic state from the interferometer. The interferometer may further comprise a plurality of waveguides arranged to pass through the interferometer to connect the M input ports to the M output ports. The plurality of waveguides may be arranged to provide a plurality of coupling locations between pairs of the plurality of waveguides. A reconfigurable beam splitter may be arranged at each of the coupling locations such that at each coupling location the two modes of electromagnetic radiation carried by the two respective waveguides are capable of coupling with each other with a reconfigurable reflection coefficient.

[0024] A reconfigurable beam splitter is understood to mean any tuneable element or device or tuneable collection of elements or devices capable of coupling two modes of electromagnetic radiation with each other with a reconfigurable effective reflection/transmission coefficient. A reconfigurable beam splitter may be capable of coupling modes with both a reconfigurable reflection coefficient and a reconfigurable phase coefficient. In some examples, a reconfigurable beam splitter may comprise a Mach-Zehnder interferometer. For example, a Mach-Zehnder interferometer may comprise two 50:50 beam splitters with a reconfigurable thermo-optic phase shifter on one internal path between the two 50:50 beam splitters. Optionally, a further phase shifter may be arranged on an external path (input or output) of a Mach-Zehnder interferometer.

[0025] The plurality of coupling locations may be arranged such that at least one of the M input modes couples with each of the other M — 1 modes in the interferometer. The interferometer may comprise fewer than — 1) /2 coupling locations.

[0026] The interferometer may be comprised in an integrated photonic circuit. The integrated circuit may further comprise one or more single photon sources. The integrated circuit may further comprise one or more detectors.

[0027] In some examples, the one or more processors may be further configured to control the number of photons in the input photonic state.

[0028] Using the boson sampler to determine a gradient of the cost function with respect to a selected parameter may comprise a) identifying first and second adjusted sets of parameter values, wherein the first adjusted set comprises the set of parameter values with the parameter value of that selected parameter adjusted positively by an amount, and wherein the second adjusted set comprises the set of parameter values with the parameter value of that selected parameter adjusted negatively by the same amount. Using the boson sampler to determine a gradient of the cost function with respect to a selected parameter may further comprise b) for each of the first and second adjusted sets: causing the boson sampler to be operated with the adjusted set of parameter values; receiving a response from the boson sampler, the response representative of an empirical probability distribution of measurement outcomes; mapping each distinct measurement outcome to a binary sequence of a plurality of binary sequences, wherein each element of a binary sequence corresponds to a measured output mode of the output photonic; determining, from the empirical probability distribution of measurement outcomes, a weighting for each binary sequence to which a measurement outcome is mapped; and evaluating the cost function for the adjusted set using the plurality of binary sequences and the weightings of the binary sequences to obtain an adjusted cost value. Using the boson sampler to determine a gradient of the cost function with respect to a selected parameter may comprise c) determining, from the adjusted cost values for the first and second adjusted sets, a gradient of the cost function with respect to that selected parameter.

[0029] Evaluating the objective function to identify a solution to the binary optimization problem may comprise comparing a function value of the evaluated objective function to a function value stored in memory. Determining a gradient of a cost function with respect to a selected parameter may include evaluating the objective function using each a set of binary sequences to provide a plurality of function values, and optionally storing in memory a binary sequence corresponding to the minimum function value found when determining a gradient. Advantageously, taking such an action helps to retain the optimal value if found while evaluating the objective function before the stopping condition is met.

[0030] The system may further comprise a second boson sampler. Advantageously, the system may then perform the hybrid quantum-classical processes described herein using a first boson sampler and a first mapping, and in parallel using a second boson sampler and a second mapping complementary to the first mapping. That is, the system may be capable of parallel processing. By taking advantage of parallel processing in this way, the system may advantageously determine a solution efficiently.

[0031] The set of single-photon sources may be configured to generate an input photonic state comprising M input modes, each input mode comprising a photon. Each binary sequence of the plurality of binary sequences may have a number L of elements that is fewer than or equal to the number of output modes M of the output photonic state.

[0032] The binary optimization problem may comprise a quadratic unconstrained binary optimization, QUBO, problem. For example, the QUBO problem may comprise a travelling salesman type problem. The binary optimization problem may comprise a portfolio optimization problem. The binary optimization problem may comprise a knapsack problem.

[0033] Some embodiments relate to a method for determining a solution to a binary optimization problem using a boson sampler. The boson sampler is operable to prepare an input photonic state comprising a plurality of input modes, generate an output photonic state comprising a plurality of output modes by interfering photons of the input photonic state, and producing measurement outcomes indicative of the presence or absence of photons in output modes of the output photonic state. The method comprises (i) initializing a set of parameter values, the set of parameter values for defining the configuration of the interferometer. The method further comprises (ii) for at least one selected parameter, using the boson sampler to determine a gradient with respect to that selected parameter. The method further comprises (iii) using the at least one determined gradient, updating the set of parameter values. The method further comprises (iv) repeating (ii) and (iii) until a stopping condition is satisfied. The method further comprises (v) after the stopping condition is satisfied, controlling the boson sampler to produce a plurality of measurement outcomes. The method further comprises (vi) mapping each distinct measurement outcome to a binary sequence of a plurality of binary sequences, wherein each element of a binary sequence corresponds to a measured output mode of the output photonic state. The method further comprises (vii) evaluating an objective function using binary sequences to which distinct measurement outcomes are mapped to identify a solution to the binary optimization problem.

[0034] Some embodiments relate to a computer-readable medium. The computer-readable medium has instructions stored thereon that, when executed by one or more processors in communication with a boson sampler, cause a method for determining a solution to a binary optimization problem to be executed.

[0035] A computer program and/or the code/instructions for performing such methods as described herein may be provided to an apparatus, such as a computer, on a computer -readable medium or computer program product. The computer -readable medium may comprise a non- transitory computer-readable medium. The computer-readable medium could be, for example, an electronic, magnetic, optical, infrared, electromagnetic, or semiconductor system, or a propagation medium for data transmission, for example for downloading the code over the Internet. Alternatively, the computer-readable medium could take the form of a physical computer-readable medium such as a semiconductor or solid-state memory, magnetic tape, a removable computer diskette, a random access memory (RAM) a read-only memory (ROM), a rigid magnetic disc, and an optical disk, such as a CD-ROM, CD-RW or DVD.

[0036] Some embodiments relate to a computing apparatus / computing device for determining a solution to a binary optimization problem. The computing apparatus comprises one or more processors. The computing apparatus further comprises one or more computer- readable storage media having stored thereon instructions that, when executed by the one or more processors while the computing apparatus is in communication with a boson sampler, cause the computing apparatus to: (i) initialize a set of parameter values, the set of parameter values for defining the configuration of the interferometer; (ii) for at least one selected parameter, use the boson sampler to determine a gradient with respect to that selected parameter; (iii) using the at least one determined gradient, update the set of parameter values; (iv) repeat (ii) and (iii) until a stopping condition is satisfied; (v) after the stopping condition is satisfied, control the boson sampler to produce a plurality of measurement outcomes; (vi) map each distinct measurement outcome to a binary sequence of a plurality of binary sequences, wherein each element of a binary sequence corresponds to a measured output mode of the output photonic state; and (vii) evaluate an objective function using binary sequences to which distinct measurement outcomes are mapped to identify a solution to the binary optimization problem.

[0037] Many modifications and other embodiments set out herein will come to mind to a person skilled in the art in light of the teachings presented herein. Therefore, it will be understood that the disclosure herein is not to be limited to the specific embodiments disclosed herein. Moreover, although the description provided herein provides example embodiments in the context of certain combinations of elements, steps and/or functions may be provided by alternative embodiments without departing from the scope of the disclosure.

Brief Description of the Drawings

[0038] Embodiments of the disclosure will now be described by way of example only, with reference to the accompanying figures, in which:

FIG. 1 shows an illustration of a quantum-classical hybrid algorithm for determining a solution to a binary optimization problem according to an example;

FIG. 2 shows a block diagram of a heterogeneous computing system for solving a binary optimization problem according to an example;

FIG. 3 shows a computing network including a heterogeneous computing system according to an example;

FIG. 4A shows an illustration of a first example of a boson sampler that may be comprised in a heterogeneous computing system;

FIG. 4B shows a second example of a boson sampler that may be comprised in a heterogeneous computing system;

FIG. 4C shows a third example of a boson sampler that may be comprised in a heterogeneous computing system; FIG. 4D shows a fourth example of a boson sampler that may be comprised in a heterogeneous computing system;

FIG. 4E shows a fifth example of a boson sampler that may be comprised in a heterogeneous computing system;

FIG. 4F shows a sixth example of a boson sampler that may be comprised in a heterogeneous computing system;

FIG. 4G shows example energy levels of an atomic system;

FIG. 4H shows an example quantum memory device;

FIG. 5 shows an example of a heterogeneous computing system with two boson samplers;

FIG. 6 shows a flow chart of an example method for performance by a heterogeneous computing system;

FIG. 7 shows a flow chart of an example method for determining a gradient of a cost function with respect to a parameter;

FIG. 8 shows an example of a heterogeneous computing system with four boson samplers;

FIG. 9A shows a graph illustrating how a cost function changed over sixty epochs when solving a first example problem;

FIG. 9B shows a graph illustrating how a cost function changed over sixty epochs when solving a second example problem;

FIG. 10A shows a graph illustrating how a cost function changed over fifty epochs when solving a third example problem;

FIG. 10B shows a histogram of the weightings of the five binary sequences that had non-zero weightings after the fifty epochs when solving the example problem of FIG. 10A; and

FIG. 11 shows a graph illustrating how a cost function changed over a number of epochs when solving a fourth example problem.

[0039] Throughout the description and the drawings, like reference numerals refer to like parts.

Detailed Description

[0040] Whilst various embodiments are described below, the disclosure is not limited to these embodiments, and variations of these embodiments may well fall within the scope of the disclosure.

[0041] The term “unitary transformation” as used herein is understood to mean substantially unitary. For example, loss or other errors may mean that a transformation imparted by a boson sampler is not an exact unitary transformation, but the skilled person will appreciate that the methods and systems described herein are also applicable in such circumstances. [0042] Many binary optimization problems are NP hard and so cannot be reliably solved on a classical computer alone. As the number L of variables increases, the number of candidate solutions increases exponentially as 2 ^L . An exhaustive search through all candidate solutions to find an optimal solution may be impractical for large numbers of binary variables.

[0043] Quantum computers / quantum information processors promise to provide solutions to many problems that are intractable on classical computers. To date, many different physical implementations and models for quantum computing have been proposed. Physical implementations proposed for quantum computing include nuclear magnetic resonance, ion traps, superconductors, quantum dots, electrons on liquid helium, solid state spectroscopy, and photonic quantum computing to name but a few. While quantum computing has advanced to the point at which “quantum supremacy” has been achieved, the technology is still in its infancy and a practical universal quantum computer - a quantum computer capable of efficiently simulating any other quantum computer - is not expected to be implemented in the near term. [0044] In part to address the lack of any powerful universal quantum computing resources, several quantum-classical hybrid algorithms have been proposed. In a quantum-classical hybrid algorithm, both classical and quantum resources are used to perform a computational task. Variational Quantum Eigensolvers (VQEs) are one such category of hybrid quantum- classical algorithms designed to compute ground state energies of quantum Hamiltonians, and have been used to find, for example, ground state energies of interacting fermionic Hamiltonians. However, VQEs still typically require quantum resources that are beyond those available currently or in the short-to-medium term. In particular, VQEs often require quantum computing resources that are able to perform a universal set of quantum logic operations on qubits, and the ability to measure quantum states in a plurality of complementary bases as required.

[0045] Described herein are methods and systems that utilise boson sampling as part of a quantum-classical hybrid algorithm to address binary optimization problems. A boson sampler is a non-universal quantum computer that relies on the interference of identical photons to generate its output. More particularly, a boson sampler comprises a network of optical components or elements (an interferometer) in which identical photons interfere with one another. As photons are quantum objects, the output of this optical network is described by a large superposition of all the possible outcomes. When a measurement is performed at the output of this network using one or more photon detectors, a single measurement outcome is realised from this superposition. If photon number resolving (PNR) detectors are used, then each measurement outcome may be described by an array or string or sequence of integers indicating how many photons were found in each output mode of the output state. If threshold detectors are used, then each measurement outcome may be described by an array or string or sequence of integers, for example a binary sequence, indicating whether photons were present or absent in each output mode of the output state. For a given input state, one can with repeated operations of the boson sampler build up an empirical probability distribution of bosonic configurations of the output state. For an introduction to boson samplers, the reader is directed to B. Gard et al. “An Introduction to Boson Sampling”, Quantum Physics ArXiv 1406.6767, and to S. Aaronson and A. Arkhipov, “The computational complexity of linear optics”, Proceedings of the forty-third annual ACM symposium on theory of computing, pages 333- 342, ACM, 2011.

[0046] To understand boson sampling, one may find it useful to think of a classical Galton board. When using a traditional Galton board, a plurality of substantially identical balls are dropped from an entrance hole one by one onto a vertical lattice of pegs, each of which randomly scatters the incoming balls onto pegs below, and the balls are collected in receptacles at the bottom of the board. Traditionally, the arrangement of the pegs on the Galton board ensures that with enough balls, the distribution of balls in each receptacle approaches the binomial distribution. If one reconfigures the arrangement of pegs, one can see a different distribution of the balls in the receptacles. If one adapts the Galton board to provide further entrance holes, then multiple balls may be dropped from an appropriate choice of the entrance holes to provide a selected input configuration. The input configuration of balls will then fall through the peg arrangement to the receptacles.

[0047] Somewhat analogously, in a boson sampler, photons may be initially arranged in some input configuration in which each photon occupies one of several input modes (analogous to multiple entrance holes on the adapted Galton board, except that in the adapted Galton board only a single ball can be dropped from any particular entrance hole at a time, while an input mode of a boson sampler may take more than one photon). The input modes may comprise, for example, spatial modes such as input ports to an interferometer (analogous to the spatial arrangement of the multiple entrance holes of the adapted Galton board), or may comprise other modes, for example temporal modes, that describe some other quantum degree of freedom of the bosonic system.

[0048] The input multimodal photonic state generated for the boson sampler is typically a product state (i.e. no quantum entanglement between input modes) comprising N photons distributed across M input modes. For illustrative purposes only, the present discussion is restricted to the case where the number of photons N is less than or equal to the number of input modes M, and where one photon is provided in any input mode. The skilled person will appreciate that the methods and systems described herein are also applicable when one or more input modes comprise more than one photon. Accordingly, in the present discussion the input state can without loss of generality be expressed as where d^' is the bosonic creation operator in the fc ^th mode.

[0049] In contrast to the balls of the Galton board, photons are quantum objects and are governed by bosonic statistics. This is starkly exemplified by the famous Hong -Ou-Mandel effect: if two single photons are incident on a 50/50 beam splitter at the same time, one in each of two input modes (input paths) of the beam splitter, then due to interference the output modes (output paths) taken by the photons become entangled. The output state can be expressed as a superposition of two configurations - one in which both photons are deflected to a first output mode and one in which both photons are deflected to the second output mode. Due to interference, other conceivable possibilities such as each output path carrying a single photon are cancelled out. The probability of the two photons being found in a particular one of the output modes is 50%.

[0050] The interferometer of a boson sampler applies some transformation to the input multimodal photonic state to produce an output multimodal photonic state. The actions of an interferometer may be modelled as a unitary transformation, although loss and other factors may mean that the transformation is not exactly unitary. In other words, the interferometer may provide a substantially unitary transformation. For example, photons passing through a beam splitter may be lost due to some kind of absorption at the beam splitter or scattering due to some mode not measured by the device. In this example, photon loss may mean that the number of photons in the input state and output state is not conserved, which means that the transformation is not exactly unitary. Thus, the term “unitary transformation” as used herein may include a transformation that is substantially unitary. Furthermore, a skilled person will appreciate that the methods and systems described herein are still applicable in circumstances where the transformation is not exactly unitary. The unitary transformation acts on the creation operators as where U (0) is a unitary matrix. The interferometer is defined by a set of parameters 0. One or more of the parameters 0 may characterise a single mode operation. For example, a parameter may characterise the phase shift imparted by a phase shifter of a passive linear photonic network. One or more of the parameters 0 may characterise a multimodal operation. For example, a parameter may characterise a transmission (or equivalently, a reflection) coefficient of a reconfigurable beam splitter in a passive linear photonic network. If the values {0} of one or more of the parameters 0 may be reconfigured then the boson sampler is said to be a reconfigurable boson sampler.

[0051] The unitary mapping transforms the input state into an output state that may be expressed as a superposition of the different possible configurations of the photons in the output modes as where C is a configuration, n ^C ' ^> is the number of photons in the j ^th output mode in configuration C, and a _c is the probability amplitude associated with configuration C. By tuning the parameter values {0}, the probability amplitudes associated with each configuration may be changed. Accordingly, a measurement using one or more photon number resolving detectors of the number of photons in each output mode can yield a measurement outcome representable as a string of integers corresponding to a configuration C. By operating the boson sampler a plurality N _s of times, it is possible to establish an empirical probability distribution of the bosonic configurations of the output state. One can expect that with many samples, the probability p _c of obtaining a measurement outcome corresponding to configuration C is approximately given by p _c = | a _c | ² .

[0052] The use of a boson sampler to address binary optimization problems may seem counterintuitive because measurement outcomes from boson samplers normally describe the number of photons in each output mode of a boson sampler, which may not be a binary value. Bosons are not qubits and so one cannot encode an arbitrary qubit Hamiltonian like in other variational quantum-classical hybrid algorithms. Furthermore, the mechanism behind creating entangled states in a linear interferometer is quite different from the mechanism behind creating entangled states in qubit circuits because entanglement between modes in a linear interferometer is very much related to boson statistics. The output photonic state produced in a boson sampler is not typically considered a multiqubit state. The same holds for the Fock measurement (a measurement of the number of photons in a mode) - it is quite dissimilar from a measurement in the qubit basis. For example, the user has little freedom as to the basis in which a measurement may be performed. Furthermore, boson sampling does not involve active error correction on the quantum device.

[0053] A quantum-classical hybrid process for identifying a binary sequence will now be described at a high level. An illustration of the process is provided in FIG. 1. The binary optimization task is to identify a binary sequence b of length L that optimizes an objective function F(bi). The length of the binary sequence / number of binary variables L is less than or equal to the number M of input modes of an input photonic state generated in a boson sampler. Measurement outcomes are produced by performing measurements on each of at least L of the M output modes of the output photonic state.

[0054] Consider, as an example, that a binary sequence of length L = M is sought and that a measurement outcome from a boson sampler indicates the number of photons in each of the M output modes. For example, the boson sampler may comprise one or more photon number resolving detectors capable of measuring the number of photons in each output mode. Each measurement outcome in such circumstances may comprise a string of integers, each integer indicating the number of bosons found in a corresponding output mode. If one samples from the boson sampler N _s times, then one may construct an empirical probability distribution of at most N _s integer strings (measurement outcomes). That is, one may determine from the plurality of measurement outcomes that the empirical probability of measuring the output state (EQ. 3) One ^can expect that as the number of measurements taken, N _s , grows large, the probability p _c of obtaining a particular integer string tends towards the square of the absolute value of the probability amplitude of the corresponding bosonic configuration | a _c | ² .

[0055] The measurement outcomes may be mapped to a binary sequences in one of several different ways.

[0056] As a first example of a suitable mapping, measurement outcomes may be mapped to binary sequences of length L = M by evaluating whether or not each integer of an integer string is greater than zero (photons present) or is equal to zero (photons absent). There are two approaches to this. According to one approach, if photons were present in the measured mode, then the corresponding element of the binary sequence has a first value, while if photons were absent from the measured mode, then the corresponding element of the binary sequence has a second value. According to the alternative approach, if photons were present in the measured mode, then the corresponding element of the binary sequence has the second value, while if photons were absent from the measured mode, then the corresponding element of the binary sequence has the first value. Taking the first value and the second value to be 1 and 0 respectively, these two mappings can be described by:

1 if bosons detected

Po- ^° ^{) =} { (EQ.4) 0 if bosons not detected 0 if bosons detected i ^: (EQ. 5) 1 if bosons not detected

[0057] According to the mapping of (EQ. 4), if the number of bosons in a particular output mode is greater than zero, then the corresponding binary digit is a one, while if the number of bosons in that particular output mode is equal to zero, then the corresponding binary digit is zero. According to the mapping of (EQ. 5), if the number of bosons in a particular output mode is greater than zero, then the corresponding binary digit is a zero, while if the number of bosons in that particular output mode is equal to zero, then the corresponding binary digit is one. Using either mapping, each measurement outcome is mapped to a binary sequence b = (b^, b^ ••• , bff^ of a plurality of binary sequences (j is 0 or 1) wherein each element of a binary sequence corresponds to a measured output mode of the output photonic state and has a value based on the presence or absence of photons in that corresponding measured output mode. Multiple measurement outcomes may be mapped to the same binary sequence. For example, the integer strings (3,0, 1,0) and (2, 0,2,0) will lead to the same binary sequence.

[0058] According to a second example of a suitable mapping, measurement outcomes may be mapped to binary sequences of length L = M by evaluating the parity of each integer of an integer string. There are two approaches to this. According to one approach, if the integer is even, then the corresponding element of the binary sequence has a first value, while if the integer is odd, then the corresponding character of the binary sequence has a second value. According to the alternative approach, if the integer is odd, then the corresponding element of the binary sequence has the first value, while if the integer is even then the corresponding character of the binary sequence has the second value. Taking the first value and the second value to be 0 and 1 respectively, these two parity mappings can be described by: p ₊ mod[n^ ^c \ 2] (EQ. 6) _: b^ = mod (EQ. 7) where © denotes binary addition.

[0059] According to the mapping of (EQ. 6), if the number of photons in a particular output mode is even, then the corresponding binary digit is a zero, while if the number of photons in that particular output mode is odd, then the corresponding binary digit is one. According to (EQ. 7), if the number of photons in a particular output mode is even, then the corresponding binary digit is a one, while if the number of photons in that particular mode is odd, then the corresponding binary digit is a zero. Using either mapping, each measurement outcome (n^\ n^\ ••• , is mapped to a binary sequence b = of a plurality of binary sequences (j is plus (+) or minus (-)) wherein each element of a binary sequence corresponds to a measured output mode of the output photonic state and has a value based on the parity of the number of photons in that corresponding measured output mode. As with the mappings of (EQ. 4) and (EQ. 5), multiple measurement outcomes may be mapped to the same binary sequence. For example, the integer strings (3,0, 1,0) and (1,2, 1,0) will lead to the same binary sequence.

[0060] After mapping, each binary sequence b has a corresponding weighting/probability that corresponds to the sum of the empirical probabilities of all measurement outcomes that map to that binary sequence:

[0061] If the size of a candidate binary sequence L is smaller than the number of input modes M, then one may optionally measure only L of the output modes of the photonic state or, alternatively, one may perform measurements on all M output modes and discard one or more elements of the corresponding measurement outcome. [0062] A cost function / penalty function, denoted E({0}), is defined based on the binary sequences to which one or more measurement outcomes map and the corresponding weightings of those binary sequences: where {0} represent a set of parameter values defining the parametrised transformation imparted by the interferometer of the boson sampler, b represents a binary sequence mapped from at least one of the measurement outcomes obtained from the boson sampler when configured with that set of parameter values, and /3 _b represents the weighting given to that binary sequence based on the (empirically derived) probabilities of all measurement outcomes that map to that that binary sequence. The terms “cost function” and “penalty function” are used herein irrespective of whether one seeks to minimize or maximize a related objective function.

[0063] The mapping of measurement outcomes to binary sequences can be performed by the classical computing resource used for implementing the hybrid process. Similarly, the weightings /3 _b of the binary sequences can also be computed using the classical resource. Furthermore, the cost function (EQ. 9) can also be evaluated using the classical resource.

[0064] In order to home in on an optimal binary sequence for the binary optimization problem, the cost function (EQ. 9) is reduced (e.g., minimized) or increased (e.g., maximized) with respect to the parameter values defining the interferometer. Any suitable method may be used for this, for example gradient descent or stochastic gradient descent. Once a stopping condition is satisfied, the boson sampler is operated using the finalised parameter values (i.e. the parameter values after no more iterations are performed), the output photonic state is sampled and the measurement outcomes are mapped to binary sequences.

[0065] Roughly speaking, by iteratively reconfiguring the boson sampler to reduce (e.g., minimise) or increase (e.g., maximize) the cost function, the number of candidate solutions to the binary optimization problem is reduced from 2 ^L to a number more manageably computable on a classical processor (e.g., the total number of objective function evaluations may be polynomial on L instead of exponential on L). Due to the probabilistic nature of quantum mechanics, even after convergence multiple binary sequences may have a high weighting and the binary sequence that reduces (e.g., minimizes) or increases (e.g., maximizes) the objective function F(b~) may or may not have the greatest weighting. Accordingly, a binary sequence is identified at least in part by evaluating the objective function for each of at least two binary sequences that are mapped to one or more measurement outcomes of the boson sampler operating with finalised parameter values.

[0066] FIG. 2 depicts a block diagram of a heterogeneous system 200 in which illustrative embodiments may be implemented. Other architectures to that shown in FIG. 2 may be used as will be appreciated by the skilled person.

[0067] System 200 is an example of a specialised computing apparatus in which computer usable program code or instructions implementing the processes may be located. In this example, system 200 includes communications fabric 202, which provides communications between processor unit(s) 204, memory unit(s) 206, input/output unit 208, communications module 210, display 212, and boson sampler 214.

[0068] The one or more processing units / processors 204 are configured to execute instructions for software that may be loaded into the memory 206. Processor unit(s) 204 may be a set of one or more processors or may be a multi-processor core, depending on the particular implementation. Furthermore, processor unit(s) 204 may be implemented using one or more heterogeneous processor systems in which a main processor is present with secondary processors on a single chip. The one or more processing unit(s) 204 may comprise one or more central processing units (CPUs), one or more graphics processing units (GPUs) or any combination thereof.

[0069] The one or more memory unit(s) 206 may comprise any piece of hardware that is capable of storing information, such as, for example, data, program code in functional form, and/or other suitable information on a temporary basis and/or a permanent basis. The one or more memory units 206 may include, for example, a random access memory or any other suitable volatile or non-volatile storage device. The one or more memory units 206 may include a form of persistent storage, for example a hard drive, a flash memory, a rewritable optical disk, a rewritable magnetic tape, or some combination thereof. The media used for persistent storage may also be removable. For example, the one or more memory units 206 may include a removable hard drive.

[0070] Input/Output unit 208 enables the input and output of data with other devices that may be in communication with the system 200. For example, input/output unit 208 may provide a connection for user input through a keyboard, a mouse, and/or other suitable devices. The input/output unit 208 may provide outputs to, for example, a printer.

[0071] Communications module 210 enables communications with other data processing systems or devices. The communications module 210 may provide communications through the use of either or both physical and wireless communications links. For example, the communications module 210 may be configured to communicate with other data processing systems or devices via a wired local area network connection, via WiFi or over a wide area network such as the internet.

[0072] Instructions for the applications and/or programs may be located in the one or more memory units 206, which are in communication with the one or more processor units 204 through communications fabric 202. Computer-implementable instructions may be in a functional form on persistent storage in the memory unit(s) 206, and may be performed by processor unit 204. These instructions may sometimes be referred to as program code, computer usable program code, or computer -readable program code that may be read and executed by a processor in processor unit 204. The program code in the different embodiments may be embodied on different physical or tangible computer-readable media.

[0073] The boson sampler 214 comprises a state generation unit 216, a linear interferometer 218, a state detection unit 220 and a dedicated (classical) control unit 222.

[0074] The state generation unit 216 is configured to generate an input photonic state |'P _/W ) comprising a plurality of input modes. The input photonic state is a product state (in other words, there is no quantum entanglement between input modes) comprising a plurality of N photons distributed across a plurality of M input modes.

[0075] The state generation module 216 may comprise a set of one or more single-photon sources. For example, a single-photon source may comprise a non-linear photonic material (such as periodically-poled lithium niobate (PPLN) or potassium titanyl phosphate (KTP)) configured to receive a pump beam from a pump laser and to probabilistically generate pairs of entangled photons, and may further comprise a photodetector configured to detect a photon of the entangled pair, thereby heralding the presence of the other photon of the pair.

[0076] The interferometer 218 comprises optical elements arranged to interfere the modes of the input photonic state, thereby transforming the input photonic state to produce an output photonic state. The interferometer 218 is configured to receive the input photonic state, to transform the input photonic state to an output photonic state, and to output the output photonic state to the state detection unit 220. The transformation is dependent on the values {0} of a set of parameters 0.

[0077] The interferometer 218 may be designed and manufactured in any suitable and desired way e.g. depending on the modes of electromagnetic radiation to be transformed by the interferometer 218. Thus, for example, when the electromagnetic radiation has an optical or infrared wavelength (e.g. between 400nm and 700nm or between 700nm and 1600nm), the optical paths through the interferometer 218 may be implemented at least partially using optical fibres. In some examples, the interferometer 218 may be implemented in bulk optics. However, in other examples, the interferometer 218 may comprise (i.e. is designed and manufactured using) an integrated circuit. In the integrated circuit, the optical paths may be implemented with, for example, a plurality of etched waveguides and plurality of coupling locations arranged in the integrated circuit. At each coupling location, tuneable elements may be arranged (e.g. EOM phase shifters) that are configured to control the coupling interaction between the waveguides. The integrated circuit may be implemented in silicon nitride (Si3N4) or any other suitable material.

[0078] The state detection unit 220 comprises one or more photon detectors configured to measure output modes of the output photonic state to produce a measurement outcome indicating whether or not photons were present in each measured output mode.

[0079] In some examples, the photon detectors may comprise photon number resolving (PNR) detectors, capable of determining how many photons are received. For example, the detectors may comprise superconducting nanowire detectors that generate an output signal intensity proportional to the (discrete) number of photons that strike a detector. The PNR detectors may comprise transition edge sensors (TESs).

[0080] In other examples, the photon detectors may comprise threshold detectors, also known as on/off detectors. Threshold detectors are not capable of determining how many photons are received but are capable of determining the presence/absence of photons in an output mode. Threshold detectors are often capable of operating at room temperature and can be used in conjunction with the mappings of (EQ. 4) and (EQ. 5).

[0081] The controller 222 is communicatively coupled to the processor unit(s) 204, the state generation module 216, the interferometer 218 and the state detection unit 220. The controller 222 may be any suitable classical computing resource for controlling the operation of the boson sampler 214. In examples, the controller 222 is implemented in a dedicated, applicationspecific processing unit. For example, the controller 222 may comprise an application -specific integrated circuit (ASIC) or an application-specific standard product (ASSP) or another domain- specific architecture (DSA). Alternatively, the controller 222 may be implemented in adaptive computing hardware (in other words, hardware comprising configurable hardware blocks / configurable logic blocks) that has been configured to perform the required functions, for example in a configured field programmable gate array (FPGA). [0082] The controller 220 is configured to receive instructions from the processor unit(s) 204, and to send measurement outcome information back to the processor unit(s) 204. More particularly, the controller is configured to, if so directed by the one or more processor unit(s) 204, configure the interferometer 218 according to a set of parameter values {0} and thereby control the transformation of the input photonic state that is implemented by the interferometer 218. For example, the controller 222 may directly send control signals that tune the effective reflectivity / transmittance of a reconfigurable beam splitter or the phase imparted by a phase shifter. The controller 222 is further configured to, if so directed by the one or more processor unit(s) 204, generate one or more control signals to cause the one or more single-photon sources to produce photons such that the state generation unit 216 produces an input photonic state. The controller 222 may optionally be able to control which input photonic state is input into the boson sampler, for example by generating one or more control signals to control the number of photons in each input mode. For example, in a photonic boson sampler in which the state generation module comprises a plurality of single-photon sources, the controller may be configured to generate one or more control signals to cause a selected number of photons to be emitted at a particular time point.

[0083] The controller 222 is further configured to receive a response from the state detection unit 220. More particularly, the controller 222 is configured to receive measurement outcomes from the photon detectors of the state detection unit 220. For example, the measurement outcomes may comprise an electrical signal from each photon detector at which a detection event occurs. In examples wherein the photon detectors are PNR detectors, the electrical signals may further be indicative of the number of photons received.

[0084] The controller 222 is further configured to communicate the response from the state detection unit 220 to the one or more processor unit(s) 204. Measurement outcomes may be aggregated together to form a response. A response may be representative of an empirical probability distribution of the measurement outcomes. The response may be received in any form. According to some examples, the response may comprise a list of N _s measurement outcomes achieved from operating the boson sampler a plurality N _s of times and taking corresponding measurements. According to other examples, the response may comprise a list of distinct measurement outcomes and their corresponding (empirically derived) probabilities. [0085] In FIG. 2, computer-readable instructions 226 are located in a functional form on computer-readable storage medium 224 that is selectively removable and may be loaded onto or transferred to system 200 for execution by processor unit(s) 204. A computer -readable storage medium may be, for example but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, device, or any suitable combination thereof. More specific examples of the computer-readable medium include the following: a portable computer diskette, a hard disk, a random access memory (RAM), a readonly memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), a portable compact disc read-only memory (CDROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer-readable storage medium may be any tangible medium that can contain or store a program for use by or in connection with an instruction execution system, apparatus, or device.

[0086] Alternatively, computer-readable instructions 226 may be transferred to system 200 from computer-readable storage medium 224 through a communications link to communications module 210 and/or through a connection to input/output unit 208. The communications link and/or the connection may be physical or wireless.

[0087] In some illustrative embodiments, computer-implementable instructions 226 may be downloaded over a network to the memory unit(s) 206 from a remote device for use with system 200. For instance, computer-implementable instructions stored in a remote server may be downloaded over a network from the server to the system 200.

[0088] The skilled person would appreciate that the architecture described above in relation to FIG. 2 is not intended to provide limitations on the computing devices with which the methods described herein may be implemented. Instead, the skilled person would appreciate that other architectures may be applicable. For example, the system 200 may include more or fewer components. For example, the system 200 may comprise more boson samplers.

[0089] The system 200 may be implemented in any of a number of ways. In some examples, the system 200 may be provided as a number of hardware modules suitable for installation in a server or computer rack (for example a conventional 19-inch server rack). For example, the processor unit(s) 204, memory unit(s) 206, input/output unit 208, and communications module 210 may be provided in a first rack-mounted hardware module, the controller 222 may be implemented in a second rack-mounted hardware module and electronically coupled to the first hardware module, the state generation unit 216 may be implemented in a third rack-mounted hardware module electronically coupled to the controller 222, the interferometer 216 may be implemented in a fourth rack-mounted hardware module electronically coupled to the controller 222 and optical fibre-connected to the state generation unit 216, and the photon detectors of the state detection unit 220 may be provided in another hardware module electronically coupled to the controller 222 and optical fibre -connected to the interferometer module.

[0090] In other examples, the system 200 may be distributed across multiple interconnected devices coupled (at least in part) over a network such as the internet.

[0091] In FIG. 3, a computing environment 300 comprises a communication network 310, one or more classical computing apparatuses 320, and a heterogeneous computing system 200-1 that may in some examples comprise the system 200 described above in relation to FIG. 2. The computing network 300 shown in FIG. 2 is intended as an example only, as would be appreciated by the skilled person. For example, the computing network 300 may comprise more or fewer classical computing apparatuses 320.

[0092] In this example, the heterogeneous computing system 330 comprises a quantum computing resource in the form of a reconfigurable boson sampler / reconfigurable boson sampling device 214 (which may be the same as the boson sampler 214 described above in relation to FIG. 2) and a classical computing resource in the form of a dedicated control system 330 (which may comprise the one or more processor unit(s) 204, one or more memory unit(s) 206 and dedicated controller 222 of FIG. 2) . The boson sampler 214 is configured to perform a boson sampling task (as will be described further below) and the dedicated control system 330 is configured to perform classical subroutines of the quantum-classical hybrid algorithm (as will be described further below) and to control or coordinate the operation of the boson sampler 214. The skilled person will appreciate that other architectures for the heterogeneous computing system 200-1 are also suitable. For example, the heterogeneous computing system 200-1 may comprise a control system 330 and a plurality of boson samplers 214 thereby enabling parallel processing of a binary optimization task. In another example, the heterogeneous computing system 200-1 may comprise a plurality of boson samplers 214, each having a respective control system 330. Furthermore, while the control system 330 and boson sampler 214 are shown as co-located entities in FIG. 3, the skilled person would appreciate that this need not be the case - for example, the control system 330 and the boson sampler 214 may be spatially separated and communicate over a private or public communication network such as communication network 310.

[0093] The one or more classical computing apparatuses 320 and the control system 330 of the heterogeneous computing system 200-1 may communicate with each other over the communication network 310. The communication network 310 may be any known type of communication network enabling wired or wireless communication between computing apparatuses 320 and the control system 330, and could be for example a Local Area Network (LAN), a Wide Area Network (WAN), or the Internet.

[0094] A user of a computing apparatus 320 may submit a binary optimization task/problem to the heterogeneous computing system 200-1 over the communication network 310. The heterogeneous computing system 200-1 may then find a solution to the submitted binary optimization problem and return a solution to the computing apparatus 320. In particular, the control system 330 and boson sampler 214 may together execute a quantum-classical hybrid process to determine a solution to the binary optimization task by processing classical tasks on the control system 330 and quantum tasks on the boson sampler 214’ . A user may alternatively interact directly with the control system 330 in order to submit the binary optimization task. In some embodiments, a computing apparatus 320 includes at least a part of the control system 330 and the control system 330 communicates with the boson sampler 214 via the communication network 310.

[0095] In some embodiments, the computing environment 300 is a cloud computing environment where boson computing is provided as a shared service to separate users. In a first example, a cloud computing service provider operates the boson sampler 214 and allows users to use the boson sampler 214. For example, a user using a computing apparatus 320, determines how to solve a binary optimization task/problem using a boson sampler, generates control instructions (e.g., parameter values) based on the determination, and transmits the control instructions to the system 200-1. In this example, the user’s computing apparatus 320 may include a part of the control system 330 or perform one or more of the control system 330 operations. In a second example, a cloud computing service provider receives a binary optimization task/problem from a user, determines how to solve the task/problem using a boson sampler, generates control instructions based on the determination, and transmits the control instructions to the system 200-1. In this example, the computing system 200-1 may be owned or operated by the service provider or another entity.

[0096] The boson sampler 214 of FIG. 2 or FIG. 3 may comprise a spatial mode interferometer, such as the boson sampler 214a illustrated in FIG. 4A. In the boson sampler of FIG. 4A, the modes of the input photonic state are spatial modes - that is, the state is defined by the number of photons in each of a plurality of spatially distinct paths. The boson sampler 214a may be implemented, at least in part, in a photonic integrated circuit.

[0097] The state generation module 216a of FIG. 4A comprises a plurality of single -photon sources 410 configured to produce single photons. One suitable photon source technology is spontaneous parametric down-conversion (SPDC). In SPDC, a non-linear crystal is pumped with a laser and, probabilistically, entangled photons are emitted (the “signal” and the “idler”). A detector (not shown in FIG. 4A) is arranged to detect the presence of the signal which, due to the entanglement, heralds the presence of a photon in the idler mode. Other photon sources may also be used, for example solid state photon sources.

[0098] The number of single -photon sources may be greater than the number M of input modes of the input photonic state |'P _/W ) in order to account for the fact that single photons may be generated only probabilistically. The state generation module 216a of FIG. 4A comprises a multiplexer 420 to route successfully generated single photons to N input ports of the M input ports of the interferometer 218a. In the example shown in FIG. 4A, the number of single photons N is equal to the number of input ports M of the interferometer 218a.

[0099] In FIG. 4A, the interferometer 218a comprises M input ports, M output ports, and a plurality of waveguides arranged to pass through the interferometer 218a to connect the M input ports to the M output ports. The plurality of waveguides are arranged to provide a plurality of coupling locations between pairs of the plurality of waveguides. A reconfigurable beam splitter 450 is arranged at each of the coupling locations such that at each coupling location the two modes of electromagnetic radiation carried by the two respective waveguides are capable of coupling with each other with a reconfigurable reflection coefficient (transmission coefficient). The reflection (transmission) coefficient of each reconfigurable beam splitter is denoted with a theta in the figure.

[0100] A parametrised/reconfigurable beam splitter is understood to mean any tuneable element or device or tuneable collection of elements or devices capable of coupling two modes of electromagnetic radiation with each other with a reconfigurable reflection/transmission coefficient and optionally a reconfigurable phase shift coefficient (not indicated in FIG. 4A). The parametrised beam splitters may be implemented in any suitable way - for example a parametrised beam splitter may comprise a Mach-Zehnder type interferometer containing a variable phase shifter in one internal path for controlling the effective beam splitter reflection coefficient of the Mach-Zehnder interferometer. The Mach-Zehnder interferometer may further comprise an external phase shifter on one external path of the Mach-Zehnder interferometer to control the relative phases of the two modes acted upon. For example, when the interferometer 218a is implemented in an integrated circuit, a reconfigurable beamsplitter may comprise a first waveguide coupling region for coupling the electromagnetic radiation modes in each waveguide, an electro -optical phase shifting element for adjusting the phase in one of the outgoing waveguides from that first coupling region, and a second waveguide coupling region for recoupling the two electromagnetic modes output from the first waveguide coupler. For example, if the interferometer 218a is implemented in bulk optics, a reconfigurable beamsplitter may comprise two 50/50 beamsplitters and a phase shifter element arranged therebetween.

[0101] The interferometer 218a may further comprise reflective elements or mirrors and/or other passive photonic elements (not shown). Accordingly, the interferometer 218a couples the single photons received at the M input ports to the plurality of M output ports based on operations defined by a set of parameter values.

[0102] The skilled person will appreciate that other architectures for the interferometer 218a may be utilised. Of course, while in the illustration the number of input and output modes is M = 4, a spatial interferometer 218a may be provided to operate on a greater or lesser number of modes.

[0103] The state detection module 220a comprises a plurality of photon number resolving (PNR) detectors 430 each arranged to receive any photons output from a corresponding output port of the interferometer 218a. For example, the PNR detectors may comprise superconducting nano wire photon detectors. A detection controller 440 (the functionality of which may be implemented by the controller 222) receives signals from the plurality of detectors 430 and can accordingly determine a measurement outcome representative of the number of photons detected at each detector. The state detection module 220a comprises one PNR detector for each of the M output modes and accordingly the measurement outcomes are representative of the number of photons measured in all output modes of the output photonic state. The skilled person would appreciate that for a binary optimization problem in which one seeks to find a binary sequence having a length L less than the number of modes M, the measurement outcomes of some of the detectors may be discarded.

[0104] Referring in particular to the arrangement of the waveguides and coupling locations in the interferometer 218a of FIG. 4A, it is notable that there is a path through the interferometer 218a from each of the input ports to each of the output ports, such that a photon in any of the input modes of the input photonic state can appear in any output mode of the output photonic state. The interferometer 218a of FIG. 4A receives as input an input photonic state having four input modes (M = 4). The coupling locations of the waveguides at which six parametrised beam splitters 450 are located corresponds (up to single mode phase shifts) to the Reck decomposition (M. Reck et al., “Experimental realization of any discrete unitary operator”, Phys. Rev. Lett., volume 73, page 58, 1994) of any unitary transformation of four modes. That is, with reconfigurable phase shifters (single mode operations), the triangular arrangement of reconfigurable beam splitters (multimodal operations) depicted in FIG. 4A can, with suitable beam splitter coefficients and phase coefficients, be used to implement any unitary transformation of an input state comprising four spatial modes to an output state comprising four spatial modes. As there are paths through the interferometer 218a such that a photon in any input port can be coupled to any output port with a suitable choice of parameter values, the specific interferometer 218a of FIG. 4A can be described as a full-depth interferometer and the boson sampler can be described as a “full-depth” boson sampler. Furthermore, every input mode interacts with every other input mode within the interferometer 218a.

[0105] The triangular arrangement of coupling locations in the Reck decomposition may be used in universal multiport interferometers. For a universal M-port interferometer, — l)/2 multimodal operations (e.g. reconfigurable beam splitters) are required. Other full-depth arrangements are also possible, for example a rectangular decomposition.

[0106] In other examples, the boson sampler 214 may be a shallow boson sampler such as that shown in FIG. 4B. FIG. 4B illustrates a second example of a boson sampler 214b that is similar to that depicted in FIG. 4A except that the interferometer 218b of FIG. 4B has fewer optical components than the interferometer 218a of FIG. 4A. In particular, the interferometer 218b of FIG. 4B comprises five parametrised beam splitters. The arrangement of waveguides and parametrised beam splitters in the interferometer 218b of FIG. 4B is similar to the first two diagonal “slices” of the triangular Reck arrangement of interferometer 218a. The interferometer 218b has fewer coupling locations than the — 1) /2 coupling locations and reconfigurable beam splitters used for a universal interferometer. As can be seen from the figure, there is no route through the interferometer 218b for photons that enter the interferometer 218b through the bottommost input port to be routed to the topmost output port, and accordingly the interferometer 218b may be described as a “shallow” (i.e. not full-depth) interferometer. The plurality of coupling locations are arranged such that at least one of the M input modes couples with each of the other M — 1 modes in the interferometer.

[0107] FIG. 4C illustrates a third example of a boson sampler 214c that is similar to that depicted in FIG. 4B except that the interferometer 218c of FIG. 4C has fewer optical components than the interferometer 218b of FIG. 4B. In particular, the interferometer 218c comprises three coupling locations having reconfigurable beam splitters. The arrangement of waveguides and parametrised beam splitters in interferometer 218c is similar to the first slice of the triangular Reck arrangement for an arbitrary unitary operation on four modes. As can be seen from the figure, there is no route through the interferometer for photons that enter the interferometer through either of the bottom two input ports to the topmost output port, and accordingly the interferometer 218c is a shallow interferometer. The plurality of coupling locations are arranged such that at least one of the M input modes couples with each of the other M — 1 modes in the interferometer.

[0108] The skilled person will appreciate that the interferometers 218a, 218b, 218c may comprise any arrangement of optical elements in a photonic integrated circuit or in bulk or fibre optics, and that the illustrations shown in FIGs. 4A to FIG. 4C are not in any way intended to be limiting.

[0109] In other examples, the boson sampler 214 may comprise a temporal mode boson sampler such as the boson sampler 214d illustrated in FIG. 4D. In the photonic boson sampler of FIG. 4D, the modes of the input multimodal bosonic state are temporal modes. In other words, the state is defined by the number of bosons in each of a plurality of temporal modes or time bins.

[0110] The state generation module 216b of FIG. 4D comprises a single-photon source 410 configured to produce a single photon in each of a plurality of time bins, so that each photon enters the temporal interferometer 218d (also referred to as a time bin interferometer), separated from the next by a time T. AS in the boson sampler 214a of FIG. 4A, the state generation module 216b may comprise further single-photon sources 410 and a multiplexer 420 in order to reliably ensure that a single photon is generated in each time period T.

[0111] In FIG. 4D a temporal interferometer comprises a plurality of temporal mode coupling devices. Each temporal mode coupling device comprises a reconfigurable beam splitter 450 and a delay line 460. A delay line 460 is arranged to connect one input port of the reconfigurable beam splitter 450 with one output port of the reconfigurable beam splitter 450. Each delay line 460 may comprise, for example, optical fibre. Each delay line 460 may have a length CT where c is the speed of light in the fibre. In this way, photons in one temporal mode may be coupled into the delay line so as to interfere with photons in the next temporal mode on a parametrised beam splitter 450. The time -bin interferometer may comprise further optical components including further optical switches.

[0112] The parameter value (e.g. transmittance) of each parametrised beam splitter 450 can be selectively tuned by a controller 222 for each time interval. Accordingly, each temporal mode coupling device can be used to implement a diagonal slice of the Reck decomposition. For example, for four input modes, the first temporal mode coupling device can be used to implement the equivalent operations of the three reconfigurable beam splitters defined by parameters 0 ₂ and 0 ₃ shown in FIG. 4A. The second temporal mode coupling device can be used to implement the equivalent operations of the two beam splitters defined by parameters 0 ₄ and 0 ₅ shown in FIG. 4A. The third temporal mode coupling device can be used to implement the equivalent operation of the beam splitter having parameter 0 ₆ in FIG. 4A. That is, the time bin interferometer 218d shown in FIG. 4D can be used to implement a full-depth circuit for four temporal modes. In general, with M — 1 delay lines and corresponding beam splitters, a time bin interferometer can implement the full-depth Reck decomposition circuit for an M mode unitary operation.

[0113] The state detection module 220b comprises a photon number resolving detector 430 arranged to receive photons output from the interferometer 218d in each temporal mode and a detection controller 440 arranged to determine the measurement outcome. The functionality of the detection controller 440 may be implemented in controller 222.

[0114] In other examples, the boson sampler 214 may comprise a shallow temporal boson sampler such as that depicted in FIG. 4E. The boson sampler 214e depicted in FIG. 4E is similar to that shown in FIG. 4D except that the time-bin interferometer comprises a single temporal mode coupling device. As discussed above, with a single temporal mode coupling device the equivalent operations of the three beam splitters defined by parameters 0 ₂ and 0 ₃ shown in FIG. 4A can be implemented. Indeed, the first slice of the Reck decomposition for an M mode transformation can be implemented.

[0115] In some examples, the boson sampler 214 may comprise a temporal mode boson sampler that uses an atomic memory, such as that depicted in FIG. 4F. The boson sampler 214f depicted in FIG. 4F is similar to that shown in FIG. 4E except that the single temporal mode coupling device of FIG. 4F comprises a quantum memory device 470. The quantum memory device 470 may be used to controllably place, for example, a photon in the first temporal mode into a superposition of being stored in the memory and being not stored in the memory, and then to controllably release any photon (or quantum state thereof) stored in the memory at the time the second temporal mode passes through, in order to interfere temporal modes in a similar way as a reconfigurable beam splitter can be used to interfere modes. A quantum memory is further described with respect to FIGs 4G-4H. Additionally, an example of a suitable quantum memory device being used in such a way is described in international patent application number PCT/GB 2021/052447 filed on 21 September 2021 in the name of ORCA Computing Limited, the content of which is incorporated herein by reference. [0116] According to an example, a quantum memory device 470 may comprise an atomic system 484 having discrete energy states (485a, 485b, 485c), which may be termed ‘energy levels’. The atomic system may comprise a single type of neutral atoms or ions or a plurality of different types of neutral atoms and/or ions. The atomic system may be an atomic ensemble comprising a plurality of atoms, such as a gas of atoms, or ions, such as a rare earth ion locked in a host medium. Such an atomic system may use energy levels associated with electrons in inner or outer shells of atoms, ions, or crystal defects. Examples of suitable atomic ensembles include vapours containing Rubidium. Other energy levels of atomic system may be used including Rydberg - type atomic systems. Additionally, or alternatively the atomic system may comprise a single neutral atom or ion, for example single 87Rb atoms in a magneto -optical- trap (MOT). The atomic system 484 of FIG. 4G is an example of a ladder-based energy system wherein the energy states of the atomic system are energy levels of valence electrons of an atomic ensemble. The discrete energy states have a first energy state (485a), a second energy state (485b) and a third energy state (485c). The third energy state 485c has a higher energy than the second energy state 485b. The second energy state 485b has a higher energy than the first energy state 485c. In the example of FIG. 4G, the first energy state 485a of the atomic system is the 5S energy level, the second energy state 485b is the 5P energy level and the third energy state 485c is the 5D energy level, however other energy levels and atomic systems may be used.

[0117] The quantum memory device 470 may be configured to receive first electromagnetic radiation comprising at least one photon (486), the photon having a field and a first frequency; the first frequency associated with a first energy. That is, the photon of a temporal mode of the input photonic state is received as first electromagnetic radiation (486) by the quantum memory device 480. The quantum memory device 480 may be further configured to receive second electromagnetic radiation (487), the second electromagnetic radiation having a second frequency. The second frequency may correspond to a second energy such that the sum of the first energy and the second energy corresponds to a resonance between the first energy state 485a and the third energy state 485c. The second energy may be different to the energy difference between the first energy state 485a and second energy state 485b of the atomic system 484. The second energy may be different to the energy difference between the second energy state 485b and third energy state 485c of the atomic system 484.

[0118] The quantum memory device 470 may be further configured to use the second electromagnetic radiation 487 to cause the photon of the temporal mode to have at least a portion of its field stored in the atomic system 484 by one or more electrons transitioning between the first energy state 485a and third energy state 485c. The quantum memory device may be further configured to receive third electromagnetic radiation (not shown), the third electromagnetic radiation having a third frequency. The third frequency may correspond to a third energy. The third energy may be less than the energy difference between the first energy state 485a and the third energy state 485c. The third energy may be different to the energy difference between the first energy state 485a and the second energy state 485b of the atomic system 484. The third energy may be different to the energy difference between the second energy state 485b and third energy state 485c of the atomic system 484. The quantum memory device 480 may be further configured to use the third electromagnetic radiation to cause at least a part of the stored field portion to be emitted from the atomic system 484.

[0119] The quantum memory device 470 may therefore be used to store photons from the first electromagnetic (EM) radiation 486 either in a deterministic manner or in a superposition state, the superposition state of a photon having components of both being stored and not stored. This is accomplished using a two (or more) photon resonance transition between the first and third energy states (485a, 485c) of the atomic system 484 wherein, individually, the first energy and the second energy correspond to off-resonant transitions of the atomic system 484. The target photon to be stored in the system may be termed the signal photon or ‘signal’ whilst the EM radiation used to controllably store or emit the signal photon may be referred to as the control pulse, or suitably ‘read-in’ or ‘read-out’ pulses. To store at least a portion of the first EM radiation 486 in the atomic system 484, the first and second EM radiation temporally overlap in the atomic system. As well as the photon from the first EM radiation being stored, a photon from the second EM radiation is also stored via the two-photon resonance.

[0120] Neither of the photon energies of the first and second electromagnetic radiation (486, 487) correspond to the atomic transitions between the first and second energy levels (485a, 485b) or between the second and third energy levels (485b, 485c). Instead, the frequencies correspond to off resonant transitions such that the photon from the first radiation corresponds to a transition from the first level 5S to an off-resonance level 485d above the 5P level. Similarly, a photon from the second EM radiation 487 corresponds to a transition from the off- resonance level 485d to the third level 5D 485c. This example is not limited to such an off- resonance level.

[0121] The off-resonance level may be above or below the second energy state 5P. Either of the transitions between a) the first energy level 485a and the off-resonance level 485d; and, b) the off resonance level 485d and the third energy level 485c, may correspond to the frequencies of any of the first, second and third radiation provided the overall resonance condition applies, which in this example is a two-photon resonance condition. The operation to store or release a photon from the atomic system may in principle involve two or more photons each from one or more sources of EM radiation. Furthermore, in this example the frequencies of the first and second EM radiation are different, however they may be the same.

[0122] Upon receiving the third EM radiation, the atomic system is stimulated to emit the stored photon to be further stored by a further photon storage function in the apparatus. Upon emitting the photon of the first EM radiation, the photon of the second EM radiation used in the two-photon resonance is also emitted.

[0123] Typically, the second and third electromagnetic radiation are pulses of EM radiation, however continuous wave (CW) operation is also feasible. The first EM radiation is a single photon of a single temporal mode of the input multimodal bosonic state. The single photon may also be part of a train of single photons temporally separated into different time bins i.e. each of the temporal modes of the input state in turn. The frequencies of the first and third EM radiation in this example are different, however they may be the same. In this example, the frequencies of the second and third EM radiation are the same. In such a case the sum of the first energy and third energy corresponds to a resonance between the third and first energy states of the atomic system. Alternatively, the second and third EM radiation frequencies (hence energies) may be different. Where the frequencies of the second and third radiation are different, the stored photon of the first EM radiation is output at a different frequency to its first frequency when it was initially input to the atomic system.

[0124] As described above, the stored photon may be deterministically stored or non- deterministically stored. This typically depends primarily on the photon flux of the second EM radiation. The higher the number of photons in the second EM signal, the greater the probability that the photon of the first EM signal is stored. The number of photons in the second and third signal affecting the storage, transmission or release of a photon may be referred to as the photon flux (the number of photons per second) of the incident EM radiation.

[0125] In circumstances where the second EM radiation has a photon flux that may either store the first EM photon or allow the first EM photon to be transmitted through the atomic system, the photon has field components in both of these states and therefore is in a quantum superposition until the superposition state collapses, for example by being measured by a detector.

[0126] When the third EM pulse is incident upon the atomic system, after the first EM pulse, the portion of the photon field being stored by the atomic system may be emitted or may remain stored. The third EM radiation may comprise a photon flux such that the field portion stored in the atomic ensemble is in a superposition of: a) being stored by the atomic ensemble; and, b) being emitted by the atomic ensemble. Thus, similar to the second EM pulse, the third EM pulse may have a photon flux that either deterministically outputs the stored portion of the photon field or places that stored portion in a quantum superposition of being continually stored in the atomic system 484 or emitted from the atomic system 484.

[0127] FIG 4H shows an example implementation of a quantum memory device 470. The apparatus of FIG. 4H is set up using bulk optics, in particular the mirrors 492, 494, and other equipment used to generate the control pulses, however portions may be implemented using optical fibre or integrated optics. Single photons pass through a dichroic mirror 492 to enter a vapour cell 490 holding an atomic ensemble 484 with the ladder-like energy level structure of FIG. 4G. In particular, the vapour may be a rubidium vapour. A further dichroic mirror 494 is disposed the other side of the vapour cell 490 such that photons may pass through the vapour cell 490 from one mirror 492 to the other 494. A strong control field 496 enters counterpropagating through the vapour cell 490, via reflection off the dichroic mirror 494 to provide the second EM radiation referred to above in relation to FIG. 4G. The counter propagating nature of the control field 496 enables phase matching conditions to be achieved. The first and second dichroic mirrors 492, 494 are transmissive to signal source wavelengths but reflective to control pulse wavelengths.

[0128] The control field 496 is generated using a source (not shown) of pulsed electromagnetic radiation at the correct optical frequency and wavelength, which for a rubidium vapour cell as used in this example is 776nm. This source may be a mode-locked titanium sapphire laser or a CW External Cavity Diode Laser followed by a Tapered Amplifier and an intensity modulator for pulse carving, for example a Pockel’s cell and/or acousto-optic modulator. Alternatively, this could be another source of pulsed electromagnetic radiation at the correct optical frequency and wavelength. The timing of the pulse generation is optimised so that the final control field pulse 496 is overlapped temporally with the target time bin in the vapour cell 490 and is ultimately selected by control electronics (not shown in the figure). When the control field 496 is temporally overlapped in the vapour cell 490 with the target time-bin, the signal photons in the time-bin are read-in and stored.

[0129] To perform a unitary operation on the time-bins of interest, which can be selected arbitrarily, a second control pulse generated from the laser and equipment (not shown) is generated, following the first control pulse, at a time such that it will temporally overlap with a further time-bin of interest in the vapour cell 490. The control electronics sets the photon flux of the second control pulse to perform the effective beam splitter operation required, for example as requested by the user or as determined from feedback from a detector. This couples the further time-bin with the time-bin previously stored in vapour cell 490 with a strength according to the photon flux of the second control pulse (and hence emulating the reflection and transmission of a beam splitter). The output mode, resulting from the interaction of the second control pulse with the vapour cell 490 the memory 470 and is directed down the chain to either be interacted with further stored time bins in subsequent memories or directed to a detector for analysis.

[0130] The skilled person will appreciate that the boson sampler examples described above in relation to FIGs 4A-4H are illustrative examples and that a boson sampler may comprise any suitable interferometer technology. Generally speaking, the quantum-classical hybrid algorithms described herein may be implemented using an M-mode boson sampler capable of a set of M — 1 multimodal operations and therefore a set of M — 1 parameter values.

[0131] One example source of errors in photonic circuits is photon loss. Photon loss may increase exponentially with the depth of a photonic circuit - that is, the more optical elements that are used to implement the transformation U (0), the more opportunity there is for photons to be lost from the boson sampler. Passive error mitigation techniques, typically based on classical postprocessing, only postpone the inevitable loss of the quantum character of an output distribution which makes the circuit efficiently simulatable classically. Thus, it may be desirable to keep a linear interferometer “shallow” (i.e. not full-depth) in order to reduce the number of components implementing the linear interferometer and thereby reducing the loss of photons.

[0132] Advantageously, the quantum-classical hybrid methods described herein can be utilised even with shallow interferometers. In a photonic circuit, this means that a binary optimization problem can be addressed with fewer photonic elements (beam splitters, phase shifters etc.) than would be used with a full-depth circuit, and therefore fewer photons are lost leading to reduced errors in computation.

[0133] As stated further above, presence/absence mappings can be used to generate binary sequences. For example, when a measurement outcome indicates that one or more photons were detected in a particular output mode, the incident event may be mapped to a first value e.g. “1”, while if no photons are detected in that particular output mode, the incident event may be mapped to a second value e.g. “0” (see (EQ 4) above). In this way, a measurement outcome may be mapped to a binary sequence for the purposes of implementing the optimization procedure described above in relation to FIG. 1. However, by using one presence/absence mapping only, not all binary sequences may be reachable by the boson sampler - for example, if photon number is conserved through the interferometer, then at least one detector will always detect the presence of photons and so the binary sequence comprising all zeros cannot be resolved (0,0,0, ... ,0) despite being a candidate solution to the binary optimization problem. In order to obviate such difficulties, the binary optimization procedure described above in relation to FIG. 1, may be repeated using the complementary presence/absence mapping (see, for example, EQ. 5) i.e. when a measurement outcome indicates that one or more photons were detected in a particular output mode, the incident event may be mapped to the second value e.g. “0”, while if no photons are detected in that particular output mode, the incident event may be mapped to the first value e.g. “1”. By using both mappings, all 2 ^M binary sequences may be reached by the boson sampler, and hence are within the solution space searched by the binary optimization procedure described herein. With a shallow boson sampler, the fraction of the candidate binary sequences that cannot be obtained with a single mapping increases, and accordingly it is beneficial to perform the binary optimization process twice, once with a first mapping (e.g. EQ. 4) and once with a second mapping (e.g. EQ. 5) in order to better determine the globally optimal solution.

[0134] FIG. 5 shows a block diagram of an alternative heterogeneous computing system 200’ ’ . The heterogeneous computing system 200-2 is suitable for identifying a globally optimal binary sequence of length L = M equal to the number of input modes / output modes of the photonic states of the boson samplers. In FIG. 5, the heterogeneous computing system 200-2 comprises a control system 330-1 and two boson samplers 214-1 and 214-2. The first boson sampler 214-1 is substantially the same as the second boson sampler 214-2. In both boson samplers 214-1, 214-2, an input photonic state comprising M modes, each containing a single photon, is generated and a parametrised unitary transformation is applied to the input state to produce an output multimodal bosonic state that is sampled. The control system 330-1 is configured to perform a binary optimization process such as that described above in relation to FIG. 1, using the two boson samplers in parallel. When processing measurement outcomes from the first boson sampler 214-1, the control system 330-1 uses the mapping p ₀ (EQ. 4). When processing measurement outcomes from the second boson sampler 214-2, the control system 330-1 uses the mapping p ₁ (EQ. 5). Advantageously, the heterogeneous computing system 200-2 is able to compare solutions found using either boson samplers and the globally optimal solution is within the binary sequences that can be found by the heterogeneous computing system 200-2.

[0135] The skilled person will appreciate that a similar consideration applies when using the parity mappings of (EQ.6) and (EQ. 7), and this will be discussed further below in relation to FIG. 8.

[0136] FIG. 6 shows a flowchart of a method 600 for determining a solution to a binary optimization problem.

[0137] The task is to identify a binary sequence b of length L to optimize an objective function F(b) that is characteristic of the binary optimization problem. The method 600 may be performed by a hybrid computing apparatus comprising a classical computing resource and a boson sampler, wherein the classical computing resource is able to, directly or indirectly, configure and communicate with the boson sampler. For illustrative purposes only, the method 600 will be described as though performed by the heterogeneous computing system 200 of FIG. 2, although the skilled person will appreciate that the method 600 is more widely applicable. Furthermore, for illustrative purposes only the method 600 will be described with reference to the presence/absence mappings of (EQ. 4) and (EQ. 5), although the skilled person will appreciate that other mappings are also applicable.

[0138] The boson sampler 214 is configured to prepare an input photonic state comprising N single photons arranged in M input modes. For the purposes of this discussion, each of the M input modes has a single photon i.e. N = M. The boson sampler 214 is further configured to perform a parametrised unitary transformation of the input state to an output state |T' _O[/7 ’(0) having M output modes. The boson sampler 214 is further configured to perform measurements on output modes of the output photonic state to produce measurement outcomes.

[0139] The method begins at 605.

[0140] At 610, the processor unit 204 select a mapping from measurement outcomes to binary sequences. For example, the mapping p ₀ described in (EQ. 4) may be selected.

[0141] At 615 the processor unit 204 selects/initialises a set of parameter values suitable for configuring the interferometer 218 of the boson sampler 214. The parameter values may characterise single mode or multimodal operations performed by the interferometer. For example, the parameter values may control a phase shift imparted by a phase shifter, or an effective reflection coefficient of a reconfigurable beam splitter. The processor unit 204 is not required to configure the boson sampler 250 with the initial parameter values.

[0142] A penalty function (EQ. 9) may be reduced (e.g., minimised) using gradient descent. The skilled person will appreciate that other minimisation methods may be utilised, for example stochastic gradient descent. In order to perform gradient descent, at each iteration one or more parameter values 0 _k of the set of parameter values {0} is transformed to 0 _k according to:

(EQ. 10) where J] is a constant referred to as the learning rate. Updating the parameter values {0} changes the empirical probability distribution that can be obtained by the boson sampler 214, which in turn affects the weights associated with each binary sequence. As can be seen in (EQ. 10) the update to a selected parameter depends on a partial derivative of the cost function with respect to that selected parameter. At 620, for at least one selected parameter, the processor unit 204 uses the boson sampler 250 to determine a gradient of a cost function with respect to that selected parameter.

[0143] According to a first example method for using the boson sampler 214 to determine the gradient of the cost function with respect to the selected parameter, on instructions from the processor unit 204, the controller 222 configures the interferometer 218 of the boson sampler 214 to be operated with the current parameter values {0} and the controller 222 initiates the use of the boson sampler with those parameter values by generating one or more control signals to cause the single-photon source(s) of the state generation module 216 to generate an input state. The boson sampler 214 accordingly generates an input state comprising a plurality of input modes, each containing a single photon. The boson sampler 214 further measures output modes of the output photonic state. The processor unit 204 then receives from the boson sampler 214 a response representative of an empirical probability distribution of measurement outcomes. The processor unit 204 maps the measurement outcomes and the corresponding probabilities to a plurality of binary sequences and their corresponding weightings using the currently selected mapping p ₀ (EQ. 4). The processor unit 204 then evaluates the cost function £■({0: 0/ })- Upon instructions from the processor unit, the controller 222 next configures the interferometer 218 of the boson sampler 214 such that the selected parameter 0 _k is shifted by a small value e while all other parameter values are maintained, and the controller 222 initiates the use of the boson sampler 214. The processor unit(s) 204 receives a response from the boson sampler 214 representative of an empirical probability distribution of measurement outcomes, maps the measurement outcomes and the corresponding probabilities to a plurality of binary sequences and their corresponding weightings, and then evaluates the cost function £({0: 0 _k + e}). The argument {0: 0 _k + e} has been used to denote that all current parameter values have been unchanged except for the selected parameter 0 _k . The numerical approximation for the gradient of the cost function with respect to the selected parameter can be used:

(EQ. 11)

[0144] To determine the gradient according to (EQ. 11) would typically require great control over all of the parameters of the boson sampler 214 as one requires the amount e to be very small for the approximation of (EQ. 11) to hold (e.g., e/0 _fc < 0.001). Furthermore, one would typically use a large number of measurement outcomes to be able to distinguish between the two similar distributions.

[0145] A second example method for using the boson sampler 214 to determine the gradient of the cost function with respect to the selected parameter is described in relation to FIG. 7.

[0146] At 710, the processor unit 204 identifies first and second adjusted sets of parameter values. In the first adjusted set, all parameter values of the set of parameter values {0} are maintained except for the selected parameter 0 _fc , which is adjusted positively by an amount s. In the second adjusted set, all parameter values of the set of parameter values {0} are maintained except for the selected parameter 0 _fc , which is adjusted negatively by the same amount s.

[0147] At 720, the processor 204 causes the boson sampler 214 to be operated with the first adjusted set of parameter values. For example, on instructions from the processor unit 204, the controller 222 may send one or more control signals to the interferometer of the boson sampler 214 to configure the tuneable elements of the interferometer 218. The controller 222 may send one or more initiation signals to start the boson sampler. At 730, the processor unit 204 receives a response from the boson sampler, the response representative of an empirical probability distribution of measurement outcomes.

[0148] At 740, the processor unit 204 maps each distinct measurement outcome to a binary sequence of a plurality of binary sequences according to the currently selected mapping. Each element of a binary sequence corresponds to a measured output mode of the output photonic state. In this example, as the currently selected mapping is p ₀ (EQ. 4), if a detection event indicates that an output mode contained photons then the corresponding element of the binary sequence is “1” while if no photons were detected in that output mode then the corresponding element of the binary sequence is “0”.

[0149] At 750, the processor unit 204 determines, from the empirical probability distribution of measurement outcomes, a weighting (EQ. 8) for each binary sequence to which a measurement outcome is mapped.

[0150] At 760, the method comprises evaluating the cost function (EQ. 9) for the first adjusted set using the plurality of binary sequences and the weightings of the binary sequences to obtain a first adjusted cost value E({0-. 0 _k + s}).

[0151] Steps 720 to 760 are repeated for the second adjusted set to obtain a second adjusted cost value E({0 0 _k — s}). The skilled person will appreciate that the order in which the first and second adjusted cost values are obtained is not important.

[0152] Once it is determined (770) that both adjusted cost values are obtained, the method proceeds to 780. At 780, the processor unit 204 determines, from the adjusted cost values for the first and second adjusted sets, a gradient of the cost function E( 0}) with respect to the selected parameter 0 _k . In particular, the gradient is evaluated from:

[0153] Determining the gradient of the cost function E({0}) with respect to the selected parameter 0 _k in this way is beneficial as the amount s by which the selected parameter 0 _k is adjusted need not be very small. Advantageously, this means that the empirical probability distributions from which the adjusted cost values E({0-. 0 _k + s}) and E({0-. 0 _k — s}) are derived can be easier to distinguish. Furthermore, the amount s may be chosen based on the physical hardware used to implement the boson sampler. For example, in a photonic system it may be more feasible to shift the selected parameter by an amount TT/2.

[0154] Whether the gradient is computed according to the first example method described above (EQ. 11) or the second example method described above (EQ. 12), evaluating the cost function (EQ. 9) may require the objective function F to be evaluated for all distinct binary sequences that are mapped to by measurement outcomes of the boson sampler 214. The processor unit 204 optionally stores in memory 206 the minimal function value of the objective function F that is found when evaluating the objective function F on these binary sequences, and stores the corresponding binary sequence. Whenever the objective function F needs to be evaluated in determining a gradient, the resulting function value can be compared with the minimum value stored in memory and replace said value if found to be more optimal.

[0155] Returning again to FIG. 6, once the one or more gradients have been determined, at 625 the processor unit 204 updates the current set of parameter values {0}.

[0156] At 630 the processor unit 204 determines whether a stopping condition has been satisfied. Any suitable stopping condition may be utilised.

[0157] For example, determining that a stopping condition has been satisfied may comprise determining that the set of parameter values have been updated a threshold number of times i.e. that the parameter values have been updated for a threshold number of iterations/epochs. The threshold number of iterations may be selected in advance by a user of the system 200.

[0158] As another example, a stopping condition may comprise a convergence criterion. Determining that a convergence criterion has been met may comprise determining that the cost function has not changed more than a threshold amount between updates of the parameter values. As another example, determining that a convergence criterion has been met may comprise evaluating a function of the weightings of the binary sequences, such as an entropic measure of the weightings, and determining that the function is less than or greater than a predetermined threshold value.

[0159] If the stopping condition has not been satisfied, then the method returns to step 620. If the stopping condition has been satisfied, then the method proceeds to step 635.

[0160] At 635, the processor unit 204, by instructing the controller 222, causes the boson sampler 214 to be operated with the finalised set of parameter values {d m} i-e. the parameter values after the stopping condition has been satisfied. At 640 a response is received from the boson sampler 214 configured with those finalised parameter values, the response representative of an empirical probability distribution of measurement outcomes.

[0161] At 645, the processor unit 204 maps each distinct measurement outcome to a binary sequence of a plurality of binary sequences according to the current mapping. Each element of a binary sequence corresponds to a measured output mode of the output photonic state, and each element has a value based on whether or not photons were detected in that output mode.

[0162] At 650, the processor unit 204 evaluates, for each of at least two binary sequences b to which a measurement outcome is mapped, the objective function F(bi) using the binary sequence to determine a corresponding function value. [0163] The gradient descent method performed does not guarantee that, after operating the boson sampler 214 with the finalised parameters, the binary sequence with the greatest weighting is the globally optimal solution to the binary optimization problem. By iteratively tuning the parameters of the boson sampler 214, the bosonic distribution sampled by the state detection module 220 is iteratively updated, causing some binary sequences to be weighted more heavily than others. However, due to the probabilistic nature of quantum measurement, the fact that binary sequences have different multiplicities (i.e. different numbers of measurement outcomes that map to them) and the fact that the boson sampler is operated only a finite number N _s of times to produce each distribution, mean that several binary sequences are likely to have weightings that are notably greater than the others. The iterative process reduces the number of candidate binary sequences of length L from 2 ^L to a number that is manageable on a classical computing resource. Accordingly, the processor unit 204 evaluates the objective function F for at least two of the binary sequences to which a measurement outcome is mapped, but preferably evaluates the objective function for several binary sequences to which a measurement outcome is mapped. The choice of which binary sequence(s) to check may be based on any suitable criteria. For example, the processor unit 204 may evaluate the objective function F for those binary sequences whose weighting is greater than a threshold value. In some examples, the processor unit 204 may evaluate the objective function for every binary sequence that corresponds to a measurement outcome received from the boson sampler operating with the finalised parameters (all binary sequences of length L having a non-zero weighting) as for high dimensional problems (large L) this is still far fewer sequences to check than the original 2 ^L candidate solutions.

[0164] The processor unit 204 may compare the minimal found function value with a function value stored previously in memory 206 and update the function value and binary sequence stored in the memory 206 if the minimal found function value is found to be more optimal than the previously stored function value. That is, the minimal found function value may be compared with the best function value found in the process of determining a gradient (at 620). [0165] At 655, the processor unit 204 determines that only one mapping has so far been used and so selects a complementary mapping (660), in this example mapping p (EQ. 5). The method then returns to step 615 (a new set of initial parameter values for the boson sampler are selected) and progresses through to step 650 as before, using the complementary mapping whenever required. In this way, the heterogeneous computing system 200 has the possibility of finding an optimal binary sequence from the full search space of 2 ^L candidate binary sequences.

[0166] At 665, the processor unit 204 identifies a binary sequence as a solution to the binary optimization problem. The binary sequence has thus been identified based at least in part on a comparison of the function values generated by operating the boson sampler with finalised parameter values, at least once for a first mapping and at least once for a second mapping complementary to the first mapping.

[0167] At 670, the method of FIG. 6 ends.

[0168] The skilled person will appreciate that the method described above in relation to FIG. 6 may be varied in several ways. For example, with two boson samplers, a first boson sampler may be used to perform steps 605 to 650 with one mapping and the second boson sampler may be used to perform steps 605 to 650 with a complementary mapping.

[0169] While the methods described above in relation to FIG. 6 and FIG. 7 utilised the mappings of (EQ. 4) and (EQ. 5), the skilled person will appreciate that in other embodiments the parity mappings of (EQ. 6) and (EQ. 7) may be used instead. When using the parity mappings, it can be advantageous to use different input photonic states.

[0170] If the number of input modes / number of output modes is of even length (M even), then both mappings p ₊ and p_ map measurement outcomes to the same plurality of binary sequences (although, of course, for any particular measurement outcome, the mappings p ₊ and p_ map that outcome to complementary binary sequences). However, the multiplicities of each binary sequence under the two mappings are usually different.

[0171] However, if the number of input modes / number of output modes is of even length (M even), then the plurality of binary sequences to which the measurement outcomes are mapped depends on whether the number of photons in the input photonic state is even or odd. In other words, if the boson sampler 214 is operated with an input state comprising M single photons distributed among M input modes: and all output modes are measured, then the measurement outcomes map to a first plurality of binary sequences of length M if the boson sampler 214 is operated with an input state comprising M — 1 single bosons distributed among M input modes:

[0172] and all output modes are measured, then the measurement outcomes map to a second plurality of binary sequences of length M different to the first plurality of binary sequences. Accordingly, in some examples and particularly those for which the length of a binary sequence is even and equal to the number of input modes, it may be advantageous to perform the method 600 at least twice: at least once with an input state comprising M single photons distributed among M output modes (EQ. 13), and at least once with an input state comprising M — 1 single photons distributed among the M input modes (EQ. 14). In this way, a globally optimal binary sequence of even length L = M may be identified.

[0173] If the number of input modes / number of output modes is of odd length (M odd) and all output modes are measured, then the two parity mappings p ₊ and p_ map measurement outcomes to complementary pluralities of binary sequences of length M. Accordingly, by performing the method 600 of FIG. 6, using the parity mappings in place of the presence/absence mappings, then the globally optimal solution to the binary optimization problem is contained within the solution space accessible by the method. Moreover, the pluralities of binary sequences to which the measurement outcomes are mapped do not depend on whether the number of photons in the input state is even or odd. Nevertheless, as for a particular problem the globally optimal solution may have a low multiplicity, it may be beneficial to repeat the method 600 twice: once with the state of (EQ. 13) and once with the state of (EQ. 14). In this way a globally optimal binary sequence of odd length L = M may be identified.

[0174] Accordingly, to find a globally optimal sequence b of length L = M to optimize a function F(b). one should perform the method 600 twice with appropriate choices of mapping and/or input multimodal state. In some examples, the processor 204 may be configured to perform steps 615 to 650 at least four times: at least once using a first parity mapping and first choice of input state, at least once using the first parity mapping and second choice of input state, at least once using the complementary mapping and first choice of input state, and at least once using the complementary mapping and second choice of input state. In this way a globally optimal binary sequence of length M may be identified and any errors emerging due to e.g. the fact that only a finite number of measurements are taken, are mitigated. [0175] FIG. 8 shows a block diagram of an alternative heterogeneous computing system 200- 3 that exemplifies this. The heterogeneous computing system 200-3 is suitable for identifying a globally optimal binary sequence of length M equal to the number of input modes / output modes of the photonic states of the boson samplers. In FIG. 8, the binary optimization system 200-3 comprises a control system 330-2 and four boson samplers 214-1, 214-2, 214-3, 214-4. The first boson sampler 214-1 and third boson sampler 214-3 are configured to operate using an input state |'P _/W ) _M (EQ. 13) and the second boson sampler 214-2 and fourth boson sampler 214-4 are configured to operate using an input photonic state |T’ _/W ) _M-1 (EQ. 14). The control system 330-2 is configured to perform the method 600 using the four boson samplers in parallel. When processing measurement outcomes from the first boson sampler 214-1 or second boson sampler 214-2, the control system 330-2 uses the mapping p ₊ such that if a measurement outcome indicates that an output mode was found to have an even number of bosons, the value of the corresponding character of the binary sequence is 0, while if the measurement outcome indicates that an output mode was found to have an odd number of bosons, the value of the corresponding character of the binary sequence is 1. When processing measurement outcomes from the third boson sampler 214-3 or fourth boson sampler 214-4, the control system 330 uses the mapping p_ such that if a measurement outcome indicates that an output mode was found to have an even number of bosons, the value of the corresponding character of the binary sequence is 1, while if the measurement outcome indicates that an output mode was found to have an odd number of bosons, the value of the corresponding character of the binary sequence is 0. Advantageously, the heterogeneous computing system 200-3 is able to identify a binary sequence faster by operating the boson samplers 214-1 to 214-4 in parallel.

[0176] Referring again to the method 600 of FIG. 6, if the candidate solutions to the binary optimization problem have a length L that is less than the number of input / output modes M, then the method 600 need only be performed once. In such circumstances measurements of all M output modes are unnecessary and instead the boson sampler 214 may be configured to perform measurements on L of the output modes and to not perform measurements on the remaining output modes. Alternatively, the boson sampler 214 may be configured to perform measurements on all M output modes, but in mapping the resulting measurement outcomes to binary sequences one or more integers/elements of the measurement outcome may be discarded. [0177] The methods described herein will now be demonstrated with reference to some further examples.

[0178] One binary optimization problem that may be addressed using the methods described herein is Quadratic Unconstrained Binary Optimization (QUBO). QUBO is usually formulated as the problem of finding the binary sequence b that minimizes an objective function F(b) of the following form: (b) = b ^T Qb (EQ. 15) where b is a binary sequence (in vector form) of length L and Q is a real, symmetric matrix of dimension L X L.

[0179] As a first example that utilised the presence/absence mappings of (EQ.4) and (EQ. 5), a quantum-classical hybrid algorithm as described herein was simulated to find a binary sequence of L = 20 binary values to minimise a randomly generated 20 X 20 symmetric matrix. More particularly, a shallow boson sampler such as that shown in FIG. 4E was simulated, with an input photonic state comprising twenty temporal input modes (M = 20) each containing a single photon (A = 20). The simulated interferometer comprised a single temporal mode coupling device comprising a reconfigurable beam splitter and a delay line. In the simulation, the numbers of photons in all twenty of the output modes were measured and accordingly each measurement outcome comprised a string of twenty integers. A set of nineteen parameter values (corresponding to nineteen transmittance values of the reconfigurable beam splitter 450 of FIG. 4E for implementing nineteen multimodal operations between consecutive modes) was optimized to minimize a cost function E({0}) (see EQ. 9) using gradient descent with a fixed learning rate of J] = 0.1 (see EQ. 10). Gradients at each iteration were computed using the method of FIG. 7. Every time the simulated boson sampler was operated, the number of measurement outcomes was taken to be N _s = 100. All parameter values were updated in every iteration, and the stopping condition used was that the iteration has been performed 60 times (60 updates).

[0180] Gradient descent was performed for both classical presence/absence mappings (EQs. 4 and 5) in turn. FIG. 9A shows how the cost function of (EQ. 9), labelled “cost” in the figure, converged over the course of the 60 iterations. On the x-axis, the iteration number is shown, and on the y-axis the cost E({0}) is plotted. In the graph, the solid curve labelled “mapping 0” corresponds to the mapping p ₀ of (EQ. 4) and the dashed curve labelled “mapping 1” corresponds to the mapping p of (EQ. 5). As can be seen in the graph, over the course of the 60 epochs, both mappings tended towards convergence, although “mapping 0” tended to convergence around a lower cost value than the “mapping 1” curve, which indicates that the optimal binary sequence was not reachable when using “mapping 1” but was reachable when using “mapping 0”.

[0181] For both mappings, if a “better” binary sequence was found when evaluating the cost function during gradient descent, this binary sequence was stored in memory. Once the stopping condition was satisfied (in this example, 60 updates of the parameter values had been performed) the boson sampler was simulated with the corresponding finalised parameter values to produce a further N _s = 100 measurement outcomes, which were then mapped to binary sequences. The objective function was evaluated for all of these 100 measurement outcomes and a comparison was made to the “best” binary sequence stored in memory to return an optimal binary sequence. For this example, the binary optimization process was able to find a 20-variable binary sequence that reduced (e.g., minimised) the objective function (in this particular example, the minimal function value was —29.3). Separately, a brute force approach verified that the minimal function value was —29.3 and that the binary sequence identified was optimal.

[0182] As the randomly generated matrix was of size 20 X 20, the number of evaluations of the objective function F(b) required for a brute force search for the solution is two to the power of twenty, which is over one million evaluations. In contrast, for the methods described herein the number of evaluations is proportional to four multiplied by the size of the parameter value set (i.e. the number of parameters to be updated) multiplied by the number of iterations performed multiplied by the number of samples N _s taken each time. The factor of four arises because there are two mappings and because two cost function values E(6 0 _k + s) and E(0 0 _k — s) are determined when determining the gradient of the cost with respect to a parameter (if one follows the method of FIG. 7). In this example, the number of evaluations of the objective function is around 456,000.

[0183] As the size of the problem grows, the difference between the number of evaluations of the objective function F(b~) using the binary optimization methods described herein and the number of evaluations required for a brute force search increases greatly. For example, with a similar setup to that used in generating FIG. 9A, a binary optimization problem characterised by the QUBO objective function of (EQ. 15) was solved for a randomly generated, symmetric matrix Q of size 50 X 50. A brute force search would require a check of two to the power of fifty different binary sequences, which is impractical on a classical computer, while the number of evaluations of the objective function required using the method of FIG. 6 is proportional to one billion one hundred and seventy six thousand. FIG. 9B shows how the cost function behaved over the course of sixty iterations. In this example, the binary optimization procedure found a binary sequence that resulted in an objective function value of —97.7.

[0184] Further examples shall now be provided that utilise the parity mappings in place of the presence/absence mappings.

[0185] It is known that finding an optimal solution to a QUBO problem is equivalent to minimizing a classical Ising Hamiltonian: constant (EQ. 16) where s _k = 2b _k — 1 is a binary variable taking a value of +1 or —1, and where the coefficients Jij and h _k take real values. That is, finding a solution to a QUBO problem for which the binary variables take values 0 or 1 can always be rephrased as finding a binary sequence of spin values + 1 or — 1 that minimizes an Ising Hamiltonian. As an example, consider the family of Ising Hamiltonians on a Mobius graph, defined by: which has a periodic boundary condition of s _L = s ₀ . The coefficient J _a is a positive real value, and the coefficient J _b is a real value.

[0186] As an example utilising parity mappings, a quantum-classical hybrid algorithm as described herein was simulated to find a binary sequence of L = 6 binary values to minimise the Ising Hamiltonian on a Mobius graph with coefficients J _a = 2 and J _b = —1.5. More particularly, a shallow boson sampler such as that shown in FIG. 4E was simulated, with an input photonic state comprising six temporal input modes (M = 6) containing either 5 or 6 bosons (A = 5,6). In the simulation, the numbers of photons in all 6 of the output modes were measured and accordingly each measurement outcome comprised a string of 6 integers. A set of 5 parameter values (corresponding to 5 transmittance values of the reconfigurable beam splitter 450 of FIG. 4E for implementing 5 multimodal operations between consecutive modes) was optimized to minimize a cost function E({0}) (see EQ. 9) using gradient descent with a fixed learning rate of p = 0.1 (see EQ. 10). Gradients at each iteration were computed using the method of FIG. 7. Every time the simulated boson sampler was operated, the number of measurement outcomes was taken to be N _s = 80. The convergence criterion used was that the iteration has been performed 50 times (50 epochs).

[0187] As the solutions to the Hamiltonian problem comprise binary sequences having values of +1 or —1, the two classical mappings p ₊ (EQ. 6) and p_ (EQ. 7) were replaced accordingly with mappings p _A and p _B respectively: p _A . mod (EQ. 18) where p _A maps an even number of photons in an output mode to the value —1 and maps an odd number of photons in an output mode to the value + 1, and where p _B maps an odd number of photons in an output mode to the value —1 and maps an even number of photons in an output mode to the value +1.

[0188] As described further above, when using a parity mapping, better results may be obtained if one uses all four combinations of input state ((EQ. 13) and (EQ. 14)) and classical mapping ((EQ. 18 and (EQ. 19)). Gradient descent was performed for each of the four possible combinations of input state and classical mapping in turn. FIG. 10A shows how each combination converged over the course of the 50 epochs. On the x-axis in FIG. 10A the epoch is shown, and on the y-axis the cost E({0}) is plotted. In the graph, “configl” represents the combination of input state having the state of (EQ. 13) and classical mapping p _A . In the graph, “config2” represents the combination of input state having the state of (EQ. 13) and classical mapping p _B . In the graph, “config3” represents the combination of input state of (EQ. 14) and classical mapping p _A . In the graph, “config4” represents the combination of input state of (EQ. 14) and classical mapping p _B . As can be seen in the graph, over the course of the 50 epochs, “config3” and “config4” converged to the same minimal value. As described above, when the number of modes of the boson sampler M is even, both parity mappings will map to the same plurality of binary sequences, but the number of photons affects which binary sequences of length L = M are accessible. [0189] For each of the four combinations, once the convergence criterion was satisfied (in this example, 50 epochs/iterations of the parameter values had been performed) the boson sampler was simulated with the corresponding converged parameter values to produce a further N _s = 80 measurement outcomes, which were then mapped to binary sequences. The weightings of those binary sequences were also calculated from the empirical probability distribution of measurement outcomes for each of the combinations. Across all four combinations, after convergence the number of binary sequences derived from (converged) measurement outcomes had been reduced to just five possible sequences. FIG. 10B shows a histogram of the weightings (from all four combinations together) of these five sequences, and the associated objective function value for these five sequences is also shown. As can be seen from the figure, three spin sequences were identified that minimize the energy of the Hamiltonian H _mob (6, 2, -1.5). These were (1,1,1, -1, -1, -1), (-1,1, 1,1, -1, -1), (-1, -1,1, 1,1, -1), which each give a value of —8.5.

[0190] A further example that utilised the parity mapping will now be described in relation to FIG. 11. In particular, the particular problem solved was a “knapsack problem”: given a set of items, each with a weight and a value, one seeks to determine which items to include in a collection so that the total weight is less than or equal to a given limit and the total value is as large as possible. More particularly, consider a list of L objects, each labelled by an index j and having a weight Wj and value V . The knapsack problem seeks to find the binary sequence b that maximises the objective function while keeping the total weight of the selected objects lower than or equal to a maximum value Wmax.

[0191] Using the methods described herein, a solution to a knapsack problem concerning 40 objects (L = 40) each having a discrete weight of 0, 1,2,3 or 4 and a real value between 0 and 1 was found with the constraint that Wmax=50. FIG. 11 shows how, for the four combinations of input state (EQs. 13 and 14) and parity mapping (EQs. 6 and 7), solutions were converged upon. After all four runs, a solution was found that had a total weight W of exactly 50, and a total value of 19.76.

[0192] For comparison, another approach to binary optimization known as simulated annealing was used as a benchmark. For many different instances of the knapsack problem, the hybrid algorithms described herein found a better solution than could be found using simulated annealing. For example, on the problem described above, the best solution found using simulated annealing gave an inferior total value of 19.25. Accordingly, it has been demonstrated that the methods described herein can beat simulated annealing in some binary optimization problems.

[0193] Variations of the methods and systems described herein are envisaged, and a few such variations are mentioned below.

[0194] The methods and systems described herein may be used in conjunction with any suitable binary optimization problem. For example, a binary optimization problem may comprise a quadratic unconstrained binary optimization (QUBO) problem. A binary optimization problem may comprise, for example, a travelling salesman problem. A binary optimization problem may comprise a portfolio optimization problem.

[0195] While in many of the examples described above, the measurement performed on the output photonic state of the boson sampler has been described as a Fock state measurement (a measurement of the number of bosons in a mode), other measurement types may be applicable. For example, it is possible to measure the parity of the number of photons in each mode. The parity operator acts on a single mode as:

H|n) = (— l) ⁿ |n) (EQ. 22) and may be implemented using, for example, homodyne detection.

[0196] Any suitable boson sampler may be utilised in conjunction with the methods and systems described herein.

[0197] In the examples described above, the learning rate was fixed. However, this is not a requirement and an adaptive learning rate.

[0198] In the examples described above in relation to FIGs. 8A to 11, all of the parameter values were updated at every iteration. The skilled person will appreciate that this is not required - at each iteration a subset of the set of parameter values may be selected for updating. [0199] In the examples described above, the number of measurement outcomes was the same N _s every time the boson sampler is used. However, the number may be varied for different uses of the sampler.

[0200] As will be appreciated by one skilled in the art, the present disclosure may be embodied as a system, method, or computer program product. Accordingly, aspects of the present disclosure may take the form of an entirely hardware embodiment, an entirely software embodiment (including firmware, resident software, micro-code, etc.) or an embodiment combining software and hardware aspects that may all generally be referred to herein as a “circuit,” “module” or “system.” Furthermore, aspects of the present disclosure may take the form of a computer program product embodied in any one or more computer-readable medium/media having computer usable program code embodied thereon.

[0201] Aspects and embodiments of the present disclosure are described herein with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems) and computer program products according to illustrative examples. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks.

[0202] The flowchart and block diagrams in the figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present disclosure. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of code, which comprises one or more executable instructions for implementing the specified logical function(s). It should also be noted that, in some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware -based systems that perform the specified functions or acts, or combinations of special purpose hardware and computer instructions.

[0203] Each feature disclosed in this specification (including any accompanying claims, abstract or drawings), may be replaced by alternative features serving the same, equivalent or similar purpose, unless expressly stated otherwise. Thus, unless expressly stated otherwise, each feature disclosed is one example only of a generic series of equivalent or similar features. The disclosure is not restricted to the details of any foregoing embodiments. The disclosure extends to any novel one, or any novel combination, of the features disclosed in this specification (including any accompanying claims, abstract and drawings), or to any novel one, or any novel combination, of the steps of any method or process so disclosed. The claims should not be construed to cover merely the foregoing embodiments, but also any embodiments which fall within the scope of the claims.