Matlab code

28 December 2013 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 anyway…

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 or later.

If you make use of this code in your research, consider including a citation to this web page in any resulting publication. Not only is it fair to give credit when you've made use of other people's work, it is important for scientific reproducability to document any code you used to help produce your results. Also, if you found this code useful, then others probably will too! Citing this page helps others find it.

I leave it to your judgement whether you feel your results made sufficient use of this code to warrant a citation; I do not insist on it. But consider whether using this code has been as helpful to your research as the least useful paper you are citing. (Typically, this sets the bar very low!). If so, you should probably include a citation to this web page.

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

Partial trace:1 1If anyone can think of an even shorter way to code these, I'd be interested in hearing about it… \(TrX(\rho \text{ or } \psi, sys, dim)\)
Partial trace over the subsystem(s) specified by (vector) \(sys\) of a density matrix \(\rho\) or state vector \(\psi\), where the subsystem dimensions are given by vector \(dim\).
Partial transpose: \(Tx(\rho,sys,dim)\)
Partial trace of a density matrix \(\rho\) with respect to the subsystem specified by \(sys\), where the subsystem dimensions are given by vector \(dim\).
Bi/tripartite partial trace:2 2The Mathematica versions can only do two sub-systems until I get round to updating them. \(TrX23(\rho,sys,dim)\),
Bi/tripartite partial transpose: \(Tx23(\rho,sys,dim)\)
Partial trace or transpose of a bipartite or tripartite density matrix \(\rho\) with respect to subsystem \(sys\), where the subsystems have dimensions specified by vector \(dim\).
Purification: \(purification(\rho)\)

Return a minimal purification of density matrix \(\rho\).

Normalise: \(normalise(\rho \text{ or } \psi)\)
Normalise the state vector \(\psi\) or density matrix \(\rho\).
Exchange subsystems: \(sysexchange(\rho \text{ or } \psi, sys, dim)\)
Exchange the order of two subsystems (specified in vector \(sys\)) of a multi-partite density matrix \(\rho\) or state vector \(\psi\), where \(dim\) specifies the dimensions of all the subsystems.
Permute subsystems: \(syspermute(\rho \text{ or } \psi, perm, dim)\)
Permuate the order of the subsystems in a multi-partite density matrix \(\rho\) or state vector \(\psi\) according to permutation \(perm\). \(dim\) specifies the dimensions of all the subsystems.
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)\).
Schmidt decomposition: \(schmidt(\psi,dim)\)
Schmidt decomposition of a pure state \(\psi\), whose subsystem dimensions are given by the vector \(dim\).
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 = \frac14 R_{ij}\sigma_i\sigma_j\).
Spin-flip: \(flip(\rho \text{ or } \psi)\)
Implements Wooters' spin-flip operation: \(\sigma_y\otimes\sigma_y\rho^*\sigma_y\otimes\sigma_y\) for density matrices, or \(\sigma_y\otimes\sigma_y\psi^*\) for state vectors.
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.
Matrix absolute value: \(absm(M)\)
Absolute value of matrix \(M\), defined as \(|M| = \sqrt{M^\dagger{}M}\).
Integer division: \(div(a,b)\)
The integer part of \(a/b\). Not really linear algebra at all, but there you go.

Channels (CP maps)

Apply channel: \(applychan(E, \rho \text{ or } \psi, rep, dim)\)
Apply channel \(E\) to density matrix \(\rho\) or state vector \(\psi\). 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.
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.
R representation: \(chan2R(\rho)\)
Convert the Choi-Jamiolkowski state representation \(\rho\) of a qubit channel to the corresponding R representation (Pauli-basis representation).
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! It's still open for \( 0 < p < 1 \) though…

Generating Random Things

Random state vector: \(randPsi(dim)\)
Generate a random state vector of dimension \(dim\), uniformly distributed according to the Haar measure.
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?).
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?).
Random unitary: \(randU(d)\)
Generate a random \(d \times d\) unitary matrix, uniformly distributed according to the Haar measure.
Random isometry: \(randV(n,m)\)
Generate a random \(n \times m\) isometry.
Random Hermitian matrix: \(randH(dim)\)
Generate a random \(d \times d\) Hermitian matrix.
Random matrix: \(randM(dim)\)
Generate a random \(d \times d\) complex matrix.

Entanglement and Entropy

von Neumann entropy: \(VNent(\rho)\)
Von Neumann entropy of density matrix \(\rho\). Note: this is not the entropy of entanglement; for that, use \(VNent(TrX(\rho,2,[dim1,dim2]))\) or \(entE\) (below).
Rényi p-entropy: \(renyi(p,\rho)\)
Rényi p-entropy of density matrix \(\rho\).
Entropy of entanglement: \(entE(\rho,dim)\)
Entropy of entanglement of a bipartite density matrix \(\rho\) 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.
Negativity: \(negativity(\rho,dim)\)
Negativity of a bipartite density matrix \(\rho\) with subsystem dimensions given by the vector \(dim\). Defined as the absolute value of twice the sum of the negative eigenvalues of \(\rho\), or 0 if all eigenvalues are positive.
Concurrence: \(concurrence(\rho \text{ or } \psi)\)
Concurrence of a 2-qubit state vector \(\psi\) or density matrix \(\rho\), defined as \(a*d-b*c\) where \(\psi=(a,b,c,d)\), or \(max[0,\lambda_1-\lambda_2-\lambda_3-\lambda_4]\) where \(\rho flip(\rho^*)\) has eigenvalues \(\lambda\).
Maximum entangled fraction: \(maxEfrac(\rho)\)
Maximum entangled (or singlet) fraction of a 2-qubit density matrix \(\rho\), defined as the maximum over all maximally entangled states \(|\phi>\) of \(<\phi|\rho|\phi>\).

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.
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.
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.
Pack/unpack variables into/from vector: \(vpack(s,...,v,...,M,...)\), \(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.

Solve dual (LMI) form SDP: \(runLMI(c,F,G,A,b)\)

Solve the following semi-definite program, given in dual (linear matrix inequality) form:

\begin{align*} &\min_x c^T x\\ &\text{subject to }\\ &\quad \Sigma_i F^(k)_i \geq G^(k)\\ &\quad Ax = b. \end{align*}

Requires a working SeDuMi installation. (SeDuMi can be made to work under Octave, too.)

Solve primal form SDP: \(runSDP(C,A,b)\)

Solve the following semi-definite program, given in primal form:

\begin{align*} &\min_X \mathrm{trace}(CX)\\ &\text{subject to }\\ &\quad X \geq 0\\ &\quad \mathrm{tr}(A^{(i)}X) = b^{(i)}. \end{align*}

Requires a working SeDuMi installation. (SeDuMi can be made to work under Octave, too.)

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.

Interesting regions: \(maxima\), \(turn\_pt\)
Functions to detect a couple of types of interesting regions, for use with \(smartplot\).

Comments

Christina

13 October 2016 I am using your partial trace function TrX, but it doesn't work when the dimensions of some systems are one. It works when the system to be traced out is before the one dimensional system, but if it is after it doesn't work.

In particular, there is an error in the fourth line of the following: 

n = length(dim);
rdim = dim(end:-1:1);
keep = [1:n];
keep(sys) = [];
dimtrace = prod(dim(sys));
dimkeep = length(p)/dimtrace;

This is because the ones from "dim" were removed but the index of "sys" was not updated. I tried to fix this, simply by updating the "sys". However, although it seems to work with a few simple examples, when I run my code (a large complicated one with endless uses of your TrX function) it doesn't work properly sometimes.

So I thought to ask you to fix it, rather than me fixing it, just to be on the safe side.

Looking forward to hear from you, and thanks for all the functions ;)

Christina

Noora Ismail

14 April 2017 This is amazing summary that I have never seen before !!!

My prayers to You!

Thanks a lot.

Eastlife

22 August 2017 Wow, Its fantastic bro. Thanks a lot! especially for the entanglement entropy and concurrence.

Roberto

24 November 2017 Comment: the new version of randRho requires randH but the file says 'nothing'.

Question: in randRho you multiply randH by 10, which in principle has no effect. Is this for numerical stability?

sreelatha

18 December 2018 Sir, This is incredible .I wanted to write a code in matlab to read a color image in quantum bits .i was searching for information .I could nt get all information but what ever I have right now is good.

sfchien

3 June 2019 Hi

Have you developed similar codes in Python?

Thanks

Chien

Toby Cubitt

17 June 2019 @sfchien #6:

No, I haven't coded Python versions of these. I'm extremely unlikely ever to do so, as I don't use Python for quantum information calculations myself (or do much numerics at all, for that matter).

For someone with decent knowledge of Numpy/Scipy and (more importantly!) a good understanding of the quantum information theory behind them, it shouldn't be too big a task. Want to take it on?

Shaohua Pan

25 April 2020 It is a great help to us !! Thanks a lot

Goutam kr Biswas

17 July 2020 Excellently useful. Thanks to the creator.

JP

16 October 2020 This is wonderful stuff. I'll be using this in my research, porting this over to python probably. I'll let you know if I build up enough functionality that it would be useful to share!

Leave a comment

All comments are moderated. Clicking submit will open your email client and let you send your comment by email. By submitting your comment you agree to license the content under a Creative Commons Attribution-ShareAlike 4.0 International License.




Creative Commons License