Maths Code

I've collected here various functions, routines, and other bits of Matlab, Octave and Mathematica code organized by topic, that might save someone, somewhere, from re-inventing the wheel. Some of them are so simple it would probably be quicker to re-code them than find this page, but since you're here...

Comments within the code should be enough to figure out what they do and how to use them (try "help <function>" from within Matlab or Octave). No guarantee they work as advertised, but I use them myself so I do correct bugs when I come across them. The Matlab code should run under both Octave and Matlab.

All the Matlab, Octave and Mathematica code linked from this page is released under the GPL license, version 2.

Quantum Information Package

For convenience, all the Matlab/Octave functions related to quantum mechanics and quantum information theory are also available in a single bundle, as well as individually (below).

Linear algebra

If anyone can think of an even shorter way to code these, I'd be interested in hearing about it...

Partial trace: TrX(ρ or ψ, sys, dim)
Partial trace over the subsystem(s) specified by (vector) sys of a density matrix ρ or state vector ψ, where the subsystem dimensions are given by vector dim.
Matlab/Octave
Partial transpose: Tx(ρ,sys,dim)
Partial trace of a density matrix ρ with respect to the subsystem specified by sys, where the subsystem dimensions are given by vector dim.
Matlab/Octave

The Mathematica versions can only do two sub-systems until I get round to updating them.

Bi/tripartite partial trace: TrX23(ρ,sys,dim)
Bi/tripartite partial transpose: Tx23(ρ,sys,dim)
Partial trace or transpose of a bipartite or tripartite density matrix ρ with respect to subsystem sys, where the subsystems have dimensions specified by vector dim.
Matlab/Octave   Mathematica
Purification: purification(ρ)
Return a minimal purification of density matrix ρ.
Matlab/Octave
Normalise: normalise(ρ or ψ)
Normalise the state vector ψ or density matrix ρ.
Matlab/Octave
Exchange subsystems: sysexchange(ρ or ψ, sys, dim)
Exchange the order of two subsystems (specified in vector sys) of a multi-partite density matrix ρ or state vector ψ, where dim specifies the dimensions of all the subsystems.
Matlab/Octave
Permute subsystems: syspermute(ρ or ψ, perm, dim)
Permuate the order of the subsystems in a multi-partite density matrix ρ or state vector ψ according to permutation perm. dim specifies the dimensions of all the subsystems.
Matlab/Octave
Tensor product: tensor(A,B,C,...)
Tensor product of arguments (generalises built-in matlab function kron). Each argument must either be a matrix, or a cell array containing a matrix and an integer. In the latter case, the integer specifies the repeat count for the matrix, e.g. tensor(A,{B,3},C) = tensor(A,B,B,B,C).
Matlab/Octave
Schmidt decomposition: schmidt(ψ,dim)
Schmidt decomposition of a pure state ψ, whose subsystem dimensions are given by the vector dim.
Matlab/Octave
Hilbert-Schmidt representation: pauli(M), unpauli(R)
Transform a 4x4 matrix M to the Hilbert-Schmidt representation in the Pauli basis, or vice-versa,
i.e. M = ¼ Rijσiσj.
Matlab/Octave
Spin-flip: flip(ρ or ψ)
Implements Wooters' spin-flip operation: σy⊗σy ρ*σy⊗σy for density matrices, or σy⊗σyψ* for state vectors.
Matlab/Octave
Commutator: comm(A,B)
Shall I insult your intelligence more than I have done already by including this on the page? Oh alright then, if you insist. This function calculates the birthday of Julius Caesar in light-years.
Matlab/Octave
Matrix absolute value: absm(M)
Absolute value of matrix M, defined as |M| = sqrt(M†M).
Matlab/Octave
Integer division: div(a,b)
The integer part of a/b. Not really linear algebra at all, but there you go.
Matlab/Octave

Channels (CP maps)

Apply channel: applychan(E, ρ or ψ, rep, dim)
Apply channel E to density matrix ρ or state vector ψ. The representation used for E is specified by rep (choi-Jamiolkowski state, purification thereof, linear operator, Kraus decomposition, Steinspring isometry, or Pauli representation of a qubit channel). The input and output dimensions of the channel are given by the vector dim. rep and dim can often be guessed automatically, so these don't necessarily need to be supplied.
Matlab/Octave
Convert channel representation: chanconv(E,from,to,dim)
Convert between different representations of a channel E. The representation used for E is specified by from, and the desired output representation by to (see applychan, above, for the possible representations). The input and output dimensions of the channel are given by the vector dim. The from and dim arguments can often be guessed automatically, so they don't necessarily need to be supplied.
Matlab/Octave
R representation: chan2R(ρ)
Convert the Choi-Jamiolkowski state representation ρ of a qubit channel to the corresponding R representation (Pauli-basis representation).
Matlab/Octave
Minimum output Rényi entropy: chanMinRenyi(p,E,rep,dim)
Numerically estimate the minimum output Rényi p-entropy of channel E specified in representation rep and with input and output dimensions given by the vector dim. Though who cares now? The minimum output entropy conjecture for p=1 is already dead!
Matlab/Octave

Generating Random Things

Random state vector: randPsi(dim)
Generate a random state vector of dimension dim, uniformly distributed according to the Haar measure.
Matlab/Octave
Random density matrix: randRho(dim)
Generate a random density matrix of dimension dim, distributed according to some measure or other (don't worry, it's dense in the set of all density matrices, which is all you really care about, right?).
Matlab/Octave
Random channel: randChan(dimA,[dimB,[dimE]],rep)
Generate a random quantum channel with specified input, and (optionally) output and environment dimensions, distributed according some other measure or other (still dense -- were you still worrying?).
Matlab/Octave
Random unitary: randU(d)
Generate a random d × d unitary matrix, uniformly distributed according to the Haar measure.
Matlab/Octave
Random isometry: randV(n,m)
Generate a random n × m isometry.
Matlab/Octave
Random Hermitian matrix: randH(dim)
Generate a random d × d Hermitian matrix.
Matlab/Octave
Random matrix: randM(dim)
Generate a random d × d complex matrix.
Matlab/Octave

Entanglement and Entropy

von Neumann entropy: VNent(ρ)
Von Neumann entropy of density matrix ρ. Note: this is not the entropy of entanglement; for that, use VNent(TrX(ρ,2,[dim1,dim2])) or entE (below).
Matlab/Octave
Rényi p-entropy: renyi(p,ρ)
Rényi p-entropy of density matrix ρ.
Matlab/Octave
Entropy of entanglement: entE(ρ,dim)
Entropy of entanglement of a bipartite density matrix ρ with subsystem dimensions given by the vector dim. The entropy of entanglement is defined as the von-Neumann entropy of the reduced density matrix of either subsystem. For bipartite pure states, all entanglement measures are asymptotically equivalent to it.
Matlab/Octave
Negativity: negativity(ρ,dim)
Negativity of a bipartite density matrix ρ with subsystem dimensions given by the vector dim. Defined as the absolute value of twice the sum of the negative eigenvalues of ρ, or 0 if all eigenvalues are positive.
Matlab/Octave
Concurrence: concurrence(ρ or ψ)
Concurrence of a 2-qubit state vector ψ or density matrix ρ, defined as a*d-b*c where ψ=(a,b,c,d), or max[0,λ1-λ2-λ3-λ4] where ρ flip(ρ*) has eigenvalues λ.
Matlab/Octave
Maximum entangled fraction: maxEfrac(ρ)
Maximum entangled (or singlet) fraction of a 2-qubit density matrix ρ, defined as the maximum over all maximally entangled states |φ> of <φ|ρ|φ>.
Matlab/Octave

Utility

Populate some useful variables: usefulvars
Populate some useful variables for playing around with quantum information. Read the contents of the file to see what goodies you get.
Matlab/Octave
Kill tinies: killtiny(X,[tolerance])
Set any elements of M that have very small magnitude to exactly zero. M can be a matrix or scalar. Real and imaginary parts are set to zero independently if necessary. Useful for tidying up matrices before displaying them.
Matlab/Octave
Pretty-print matrices: pprintm(M,[rowlabels,[collabels]]
Pretty-print a matrix in a format that makes it easier to visually parse the structure of the matrix. Rows and columns can be labelled, and labellings suitable for state vectors or density matrices of arbitrary combinations of qudits can be generated automatically.
Matlab/Octave
Pack variables into vector: vpack(s,...,v,...,M,...)
Unpack variables from vector: vunpack(V,s,...,v,...,A,...)
Pack or unpack scalars s, vectors v, and matrices M into a single vector, or unpack them again. Arguments to vpack can be in any order. Argument V to vunpack is the vector to be unpacked, the other arguments are used to determine the dimensions of the objects to unpack into (contents of these arguments are ignored).

Intended for use with optimization routines which (used to?) require all optimization parameters to be passed in a single vector.

Matlab/Octave
Solve dual (LMI) form SDP: runLMI(c,F,G,A,b)
Solve the following semi-definite program, given in dual (linear matrix inequality) form:
minimize  cT x
subject to Σi F(k)i >= G(k)
Ax = b
Requires a working SeDuMi installation. (SeDuMi can be made to work under Octave, too.)
Matlab/Octave
Solve primal form SDP: runSDP(C,A,b)
Solve the following semi-definite program, given in primal form:

minimize  trace(CX)
subject toX >= 0
trace(A(i)X) = b(i)
Requires a working SeDuMi installation. (SeDuMi can be made to work under Octave, too.)
Matlab/Octave

Plotting

`Smart' plotting: smartplot(...)
Read the internal documentation for the low-down on this one (try "help smartplot" at a Matlab/Octave prompt). Essentially, smartplot returns points sampled from a user-supplied function, which can then be plotted with the built-in plotting functions.

The `smart' part is that a higher density of points is sampled around interesting regions (e.g. turning points). Interesting regions are identified by a user-defined function, so can be anything.

Matlab/Octave
Interesting regions: maxima, turn_pt
Functions to detect a couple of types of interesting regions, for use with smartplot.
Matlab/Octave