#####################################################################################
## MATLAB Script for Estimating the CB Norm of a Map
## Algorithm Designed by: Nathaniel Johnston, David Kribs and Vern Paulsen
## Algorithm Coded by: Nathaniel Johnston
## Script Hosted at: www.nathanieljohnston.com/cb_norm_script.asp
## Script Version: 1.00
## Last Modified: Monday, February 18, 2008
##
## Main Function: CBNorm
## Auxilliary Functions: CellMatricize, IsCPMap, MakeLinIndep, RandomPositive
#####################################################################################

with(LinearAlgebra):


#####################################################################################
## CBNorm
##
## Inputs: CelA - An array of matrices that are the left generalized Choi-Kraus
##                operators of the map
##         CelB - An array of matrices that are the right generalized Choi-Kraus
##                operators of the map
##         NumIts - The number of random positive-definite matrices to generate
##                  when trying to approximate the CB norm
##         NumOps - The number of left (or right) generalized Choi-Kraus operators
##
## Output: An upper bound for the CB norm of the given map. A higher value of NumIts
##         results in a closer (on average) estimate of the CB norm of the given map,
##         but also results in a longer computation time
#####################################################################################

CBNorm := proc(CelA,CelB,NumIts,NumOps)
    local i,x,n,k,CBGuess,NewCBGuess,ST,G,H,GG,HH,GGCB,HHCB,EF:
    n:=ColumnDimension(CelA[1]):k:=RowDimension(CelA[1]):
    CBGuess:=infinity:
    if IsCPMap(CelA,CelB,NumOps,n,k) then
        CBGuess:=Norm(add(MatrixMatrixMultiply(CelA[x],CelB[x]),x=1..NumOps),2):
    else
        EF:=MakeLinIndep(CelA,CelB,NumOps,n,k):
        for i from 1 to NumIts do
            ST:=RandomPositive(EF[3],0,1):
            for x from 1 to EF[3] do
                H[x]:=Matrix(n,k,(i,j) -> add(ST[1][x,l]*EF[2][l][i,j],l=1..EF[3])):
                G[x]:=Matrix(k,n,(i,j) -> add(ST[2][l,x]*EF[1][l][i,j],l=1..EF[3])):
                HH[x]:=MatrixMatrixMultiply(HermitianTranspose(H[x]),H[x]):
                GG[x]:=MatrixMatrixMultiply(G[x],HermitianTranspose(G[x])):
            od:
            HHCB:=simplify(Matrix(k,k,(i,j) -> add(HH[l][i,j],l=1..EF[3]))):
            GGCB:=simplify(Matrix(k,k,(i,j) -> add(GG[l][i,j],l=1..EF[3]))):
            NewCBGuess:=Re(sqrt(evalf(Norm(GGCB,2)) * evalf(Norm(HHCB,2)))):
            if NewCBGuess < CBGuess then
                CBGuess:=NewCBGuess:
            fi:
        od:
    fi:
    RETURN(CBGuess):
end:


#####################################################################################
## CellMatricize
##
## Inputs: Cel - An array of matrices
##         a - The number of matrices in the array Cel
##         b - The row dimension of the input matrices
##         c - The column dimension of the input matrices
##
## Output: An a-by-(b*c) matrix with the ith row as the ith input matrix
##         The first c entries of the ith row are the first row of the ith matrix,
##         the next c entries of the ith row are the second row of the ith matrix,
##         and so on
#####################################################################################

CellMatricize := proc(Cel,a,b,c)
    RETURN(Matrix(a,b*c,(i,j) -> Cel[i][floor((j-1)/c)+1,j-c*floor((j-1)/c)])):
end:



#####################################################################################
## IsCPMap
##
## Inputs: CelA - An array of matrices that are the left generalized Choi-Kraus
##                operators of the map
##         CelB - An array of matrices that are the right generalized Choi-Kraus
##                operators of the map
##         NumOps - The number of left (or right) generalized Choi-Kraus operators
##         n - The dimension of the matrix space that the given map maps from
##         k - The dimension of the matrix space that the given map maps to
##
## Output: Either true or false, indicating whether or not the given map is
##         completely positive
#####################################################################################

IsCPMap := proc(CelA,CelB,NumOps,n,k)
    local i,j,x,Choi,LgChoi,MtxUnit:
    if not (n = k) then
        RETURN(false):
    else
        for i from 1 to n do
            for j from 1 to n do
                MtxUnit[i][j]:=OuterProductMatrix(UnitVector(i,n),UnitVector(j,n)):
                Choi[i][j]:=add(MatrixMatrixMultiply(MatrixMatrixMultiply(A[x],MtxUnit[i][j]),B[x]),x=1..NumOps):
            od:
        od:
        LgChoi:=Matrix(n^2,n^2,(i,j) -> Choi[floor((i-1)/n)+1][floor((j-1)/n)+1][((i-1) mod n)+1,((j-1) mod n)+1]):
        RETURN(IsDefinite(LgChoi)):
    fi:
end:


#####################################################################################
## MakeLinIndep
##
## Inputs: CelA - An array of matrices that are the left generalized Choi-Kraus
##                operators of the map
##         CelB - An array of matrices that are the right generalized Choi-Kraus
##                operators of the map
##         NumOps - The number of left (or right) generalized Choi-Kraus operators
##         n - The dimension of the matrix space that the given map maps from
##         k - The dimension of the matrix space that the given map maps to
##
## Outputs: (1) - A linearly independent array of left generalized Choi-Kraus
##                operators for the input map
##          (2) - A linearly independent array of right generalized Choi-Kraus
##                operators for the input map
##          (3) - The number of linearly independent left (or right) generalized
##                Choi-Kraus operators
#####################################################################################

MakeLinIndep := proc(CelA,CelB,NumOps,n,k)
    local BM,BS,u,x,i,j,bg,BG,d,C,DOp,CM,CS,v,cg,CG,c,E,F:
    BM:=CellMatricize(CelB,NumOps,n,k):
    BS:=Basis({Row(BM,[1..NumOps])}):
    u:=round(Rank(BM)):
    for x from 1 to NumOps do
        for i from 1 to n*k do
            for j from 1 to u do
                bg[i,j] := BS[j][i]:
            od:
            bg[i,u+1] := BM[x,i]:
        od:
        BG[x]:=ReducedRowEchelonForm(Matrix(n*k,u+1,(i,j) -> bg[i,j])):
        for j from 1 to u do
            d[x,j]:=BG[x][j,u+1]:
        od:
    od:
    for x from 1 to u do
        C[x]:=Matrix(k,n,(i,j) -> add(d[l,x]*CelA[l][i,j],l=1..NumOps)):
        DOp[x]:=Matrix(n,k,(i,j) -> BS[x][j+(i-1)*k]):
    od:
    CM:=CellMatricize(C,u,k,n):
    v:=round(Rank(CM)):
    CS:=Basis({Row(CM,[1..u])}):
    for x from 1 to u do
        for i from 1 to n*k do
            for j from 1 to v do
                cg[i,j] := CS[j][i]:
            od:
            cg[i,v+1] := CM[x,i]:
        od:
        CG[x]:=ReducedRowEchelonForm(Matrix(n*k,v+1,(i,j) -> cg[i,j])):
        for j from 1 to v do
            c[x,j]:=CG[x][j,v+1]:
        od:
    od:
    for x from 1 to v do
        E[x]:=Matrix(k,n,(i,j)->CS[x][j+(i-1)*n]):
        F[x]:=Matrix(n,k,(i,j)->add(c[l,x]*DOp[l][i,j],l=1..u)):
    od:
    RETURN(E,F,v):
end:


#####################################################################################
## RandomPositive
##
## Inputs: ndim - the number of rows (or columns) for the matrix to have
##         evalLower - the smallest eigenvalue allowable for the random matrix
##         evalUpper - the largest eigenvalue allowable for the random matrix
##
## Outputs: (1) - An ndim-by-ndim positive-definite (not necessarily hermitian)
##                matrix with eigenvalues in the interval (evalLower, evalUpper]
##          (2) - The inverse matrix of (1)
#####################################################################################

RandomPositive := proc(ndim,evalLower,evalUpper)
    local r,i,RD,RM,V,W,U,P:
    r:=0:
    RD:=RandomMatrix(ndim,ndim,generator=evalLower+DBL_EPSILON..evalUpper,outputoptions=[shape=diagonal]) + RandomMatrix(ndim,ndim,generator=evalLower..evalUpper,outputoptions=[shape=skewsymmetric])+I*RandomMatrix(ndim,ndim,generator=evalLower..evalUpper,outputoptions=[shape=symmetric]):
    while r < ndim do
        RM:=RandomMatrix(ndim,ndim,generator=0.0..1.0) + I*RandomMatrix(ndim,ndim,generator=0.0..1.0):
        r:=round(Rank(RM)):
    od:
    for i from 1 to ndim do
        V[i]:=Vector(ndim,(j) -> RM[j,i]):
    od:
    W:=GramSchmidt(convert(V,list),normalized):
    U:=Matrix(ndim,ndim,(i,j) -> W[j][i]):
    P[1]:=MatrixMatrixMultiply(MatrixMatrixMultiply(U,RD),HermitianTranspose(U)):
    P[2]:=MatrixInverse(P[1]):
    RETURN(P[1],P[2]):
end:


#####################################################################################
## After loading in the previous procedures, you may use code similar to the
## following to estimate the completely bounded norm of a map.
#####################################################################################

> NumOps := 2;
> NumIts := 100;
> A := Array(1 .. NumOps);
> B := Array(1 .. NumOps);
> A[1] := DiagonalMatrix([exp(-(1/4)*(5*I)*Pi), exp(-I*Pi), exp(-(1/4)*(3*I)*Pi)]);
> A[2] := IdentityMatrix(3);
> B[1] := DiagonalMatrix([exp((1/4)*(5*I)*Pi), exp(I*Pi), exp((1/4)*(3*I)*Pi)]);
> B[2] := -IdentityMatrix(3);
> CBNorm(A, B, NumIts, NumOps);