fork download
copy 
  1. import cvxpy as cp
  2. import numpy as np
  3. import scipy.linalg as la
  4.  
  5. def fix_density_matrix_max_entropy(rho, max_iters=1000, trace_penalty=1e6, positivity_penalty=1e6):
  6.     """
  7.     Fixes the off-diagonal elements of a density matrix and maximizes its entropy.
  8.     rho: Input density matrix (numpy array)
  9.     max_iters: Maximum number of iterations for the optimization solver
  10.     trace_penalty: Penalty coefficient for the trace constraint
  11.     positivity_penalty: Penalty coefficient for the positivity constraint
  12.     """
  13.     n = rho.shape[0]
  14.  
  15.     # Define eigenvalues as cvxpy variables
  16.     eigenvalues = cp.Variable(n, nonneg=True)
  17.  
  18.     # Define the Hermitian density matrix as a cvxpy variable
  19.     rho_var = cp.Variable((n, n), complex=True)
  20.  
  21.     # Ensure the matrix is Hermitian
  22.     constraints = [rho_var == cp.conj(rho_var.T)]
  23.  
  24.     # Apply fixed off-diagonal elements
  25.     for i in range(n):
  26.         for j in range(n):
  27.             if i != j:  # Fix off-diagonal elements
  28.                 constraints.append(rho_var[i, j] == rho[i, j])
  29.  
  30.     # Constraint to ensure positive semidefiniteness (explicitly enforce positive semidefiniteness)
  31.     constraints.append(rho_var >> 0)
  32.  
  33.     # Constraint on eigenvalues
  34.     constraints.append(cp.sum(eigenvalues) == 1)  # Sum of eigenvalues must equal 1
  35.  
  36.     # Small value for numerical stability
  37.     epsilon = 1e-10
  38.  
  39.     # Penalty term for the trace constraint
  40.     trace_error = cp.abs(cp.trace(rho_var) - 1)
  41.  
  42.     # Penalty term for positive semidefiniteness
  43.     # Add a penalty if the smallest eigenvalue is negative
  44.     min_eigenvalue = cp.lambda_min(rho_var)
  45.     positivity_error = cp.maximum(-min_eigenvalue, 0)
  46.  
  47.     # Objective function: Maximize entropy + penalties for constraints
  48.     obj = cp.Minimize(
  49.         cp.sum(-cp.entr(eigenvalues + epsilon)) +
  50.         trace_penalty * trace_error +
  51.         positivity_penalty * positivity_error
  52.     )
  53.  
  54.     # Define the optimization problem
  55.     prob = cp.Problem(obj, constraints)
  56.  
  57.     # Suppress debug output by setting verbose to False
  58.     try:
  59.         prob.solve(solver=cp.SCS, verbose=False, max_iters=max_iters)
  60.     except cp.error.SolverError:
  61.         # If SCS solver fails, try MOSEK (if available)
  62.         try:
  63.             prob.solve(solver=cp.MOSEK, verbose=False)
  64.         except cp.error.SolverError:
  65.             raise ValueError("Optimization failed.")
  66.  
  67.     if prob.status not in ["optimal", "optimal_inaccurate"]:
  68.         raise ValueError(f"Optimization failed: {prob.status}")
  69.  
  70.     return rho_var.value
  71.  
  72. def generate_density_matrix(n, seed, max_eigenvalue=1.0):
  73.     """Generates a random density matrix."""
  74.     np.random.seed(seed)
  75.  
  76.     while True:
  77.         # Generate a random Hermitian matrix
  78.         a = np.random.rand(n, n) + 1j * np.random.rand(n, n)
  79.         h = (a + a.conj().T) / 2  # Define h outside the loop
  80.  
  81.         # Ensure positive definiteness using the exponential map
  82.         rho = la.expm(h)
  83.         rho /= np.trace(rho)
  84.  
  85.         # Convert diagonal elements to real numbers
  86.         np.fill_diagonal(rho, np.real(np.diag(rho)))
  87.  
  88.         # Limit the maximum eigenvalue
  89.         eigenvalues = np.linalg.eigvals(rho)
  90.         if np.max(np.abs(eigenvalues)) <= max_eigenvalue:
  91.             return rho
  92.  
  93. def main(seeds, dimensions):
  94.     """
  95.     Fixes density matrices for multiple seeds and dimensions, and outputs the results.
  96.     seeds: List of seed values
  97.     dimensions: List of dimensions
  98.     """
  99.     for n in dimensions:
  100.         for seed in seeds:
  101.             print(f"\n=== Final Result, Dimension {n}, Seed {seed} ===")
  102.  
  103.             # Generate the original density matrix
  104.             rho = generate_density_matrix(n, seed)
  105.  
  106.             # Compute the corrected density matrix
  107.             rho_fixed = fix_density_matrix_max_entropy(rho)
  108.  
  109.             # Compute eigenvalues
  110.             eigenvalues = np.linalg.eigvals(rho_fixed)
  111.  
  112.             # Check positive semidefiniteness
  113.             is_positive_semidefinite = np.all(eigenvalues >= -1e-10)
  114.  
  115.             # Check the trace
  116.             trace_fixed = np.trace(rho_fixed)
  117.  
  118.             # Maximum error
  119.             max_error = np.max(np.abs(rho_fixed - rho))
  120.  
  121.             # Output the results
  122.             print("Original Density Matrix:")
  123.             print(rho)
  124.             print("\nCorrected Density Matrix:")
  125.             print(rho_fixed)
  126.             print("\nEigenvalues:", eigenvalues)
  127.             print("Positive Semidefiniteness Check:", is_positive_semidefinite)
  128.             print(f"Trace: {trace_fixed:.15f}")
  129.             print("Maximum Error:", max_error)
  130.  
  131. # Example execution
  132. seeds = [1]  # List of seed values
  133. dimensions = [50]  # List of dimensions
  134. main(seeds, dimensions)
Runtime error #stdin #stdout #stderr 0.14s 23308KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
Traceback (most recent call last):
  File "./prog.py", line 1, in <module>
    import cvxpy as cp
ModuleNotFoundError: No module named 'cvxpy'
https://ideone.com/56hmuK
language:
Python 3 nbc (python 3.7.3)
created:
2 days ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!