Commit afa34632 authored by Guillaume Sagnol's avatar Guillaume Sagnol


parent 5bf4ce02
......@@ -33,6 +33,7 @@ import numpy as np
import sys
import six
from six.moves import zip, range
import itertools
from .tools import *
from .constraint import *
......@@ -694,6 +695,71 @@ class AffinExp(Expression):
cf. doc of :func:`picos.partial_transpose() <>`
def partial_trace(self,k=1,dim=None):
"""partial trace
cf. doc of :func:`picos.partial_trace() <>`
sz = self.size
if dim is None:
if sz[0] == sz[1] and (sz[0]**0.5) == int(sz[0]**0.5) and (sz[1]**0.5) == int(sz[1]**0.5):
dim = (int(sz[0]**0.5),int(sz[1]**0.5))
raise ValueError('The default parameter dim=None assumes X is a n**2 x n**2 matrix')
pdim = np.product(dim)
if pdim != sz[0] or pdim != sz[1]:
raise ValueError('The product of the sub-dimensions does not match the size of X')
if k > len(dim):
raise Exception('we must have k <= len(dim)')
dim_reduced = [d for d in dim]
del dim_reduced[k]
dim_reduced = tuple(dim_reduced)
pdimred = np.product(dim_reduced)
fact = cvx.spmatrix([],[],[],(pdimred**2, pdim**2))
for iii in itertools.product(*[range(i) for i in dim_reduced]):
for jjj in itertools.product(*[range(j) for j in dim_reduced]):
#element iii,jjj of the partial trace
row = int(sum([iii[j] * np.product(dim_reduced[j+1:]) for j in range(len(dim_reduced))]))
col = int(sum([jjj[j] * np.product(dim_reduced[j+1:]) for j in range(len(dim_reduced))]))
#this corresponds to the element row,col in the matrix basis
rowij = col * pdimred + row
#this corresponds to the elem rowij in vectorized form
#computes the partial trace for iii,jjj
for l in range(dim[k]):
iili = list(iii)
iili = tuple(iili)
jjlj = list(jjj)
jjlj = tuple(jjlj)
row_l = int(sum([iili[j] * np.product(dim[j+1:]) for j in range(len(dim))]))
col_l = int(sum([jjlj[j] * np.product(dim[j+1:]) for j in range(len(dim))]))
colij_l = col_l * pdim + row_l
fact[rowij,colij_l] = 1
newfacs = {}
for x in self.factors:
newfacs[x] = fact * self.factors[x]
if self.constant:
cons = fact * self.constant
cons = None
return AffinExp(newfacs,cons,(pdimred,pdimred),'Tr_'+str(k)+'('+self.string+')')
def hadamard(self,fact):
"""hadamard (elementwise) product"""
return self^fact
......@@ -72,6 +72,7 @@ __all__=['_retrieve_matrix',
......@@ -546,6 +547,45 @@ def partial_transpose(exp,dim = None):
return exp.partial_transpose(dim)
def partial_transpose(X,k=1,dim = None):
r"""Partial trace of an Affine Expression,
with respect to the ``k``th subsystem for a tensor product of dimensions ``dim``.
If ``X`` is matrix
:class:`AffinExp <picos.AffinExp>`
that can be written as :math:`X = A_0 \otimes \cdots \otimes A_{n-1}`
for some *square matrices* :math:`A_0,\ldots,A[n-1]`
of respective size ``dim[0]``, ..., ``dim[n-1]``,
this function returns the matrix
:math:`Y = \operatorname{trace}(A_k)\quad A_0 \otimes \cdots A_{k-1} \otimes A_{k+1} \otimes \cdots \otimes A_{n-1}`.
The default value ``dim=None`` automatically computes the size of the subblocks,
assuming that ``X`` is a :math:`n^2 \times n^2`-square matrix
with blocks of size :math:`n \times n`.
#TODO !!!
>>> import picos as pic
>>> import cvxopt as cvx
>>> P = pic.Problem()
>>> X = P.add_variable('X',(4,4))
>>> X.value = cvx.matrix(range(16),(4,4))
>>> print X #doctest: +NORMALIZE_WHITESPACE
[ 0.00e+00 4.00e+00 8.00e+00 1.20e+01]
[ 1.00e+00 5.00e+00 9.00e+00 1.30e+01]
[ 2.00e+00 6.00e+00 1.00e+01 1.40e+01]
[ 3.00e+00 7.00e+00 1.10e+01 1.50e+01]
>>> print pic.partial_trace(X) #partial transpose with respect to second subsystem (k=1) #doctest: +NORMALIZE_WHITESPACE
[ 5.00e+00 2.10e+01]
[ 9.00e+00 2.50e+01]
>>> print pic.partial_trace(X,0) #and now with respect to first subsystem (k=0) #doctest: +NORMALIZE_WHITESPACE
[ 1.00e+01 1.80e+01]
[ 1.20e+01 2.00e+01]
return X.partial_trace(k,dim)
def detrootn(exp):
"""returns a :class:`DetRootN_Exp <picos.DetRootN_Exp>` object representing the determinant of the
:math:`n` th-root of the symmetric matrix ``exp``, where :math:`n` is the dimension of the matrix.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment