main > main: presentation for paper-dilation
This commit is contained in:
parent
669796b805
commit
88e539929b
|
|
@ -0,0 +1,229 @@
|
|||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Examples - non-dilatable $d$-parameter $C_{0}$-semigroups #\n",
|
||||
"\n",
|
||||
"This notebook provides supplementary material to the research paper <https://doi.org/10.48550/arXiv.2210.02353>.\n",
|
||||
"\n",
|
||||
"In §5.3 a construction is given which yields for any $d \\geq 2$\n",
|
||||
"a $d$-parameter $C_{0}$-semigroup, $T$ satisfying:\n",
|
||||
"\n",
|
||||
"- $T$ is contractive (equivalently its marginals $T_{i}$ are for each $i$);\n",
|
||||
"- the generator $A_{i}$ of $T_{i}$ has strictly negative spectral bound for each $i$;\n",
|
||||
"\n",
|
||||
"and such that $T$ is __not__ **completely dissipative**.\n",
|
||||
"Thus by the classification Theorem (Thm 1.1), $T$ does not have a regular unitary dilation.\n",
|
||||
"\n",
|
||||
"This Notebook demonstrates this general result empirically."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import os;\n",
|
||||
"import sys;\n",
|
||||
"\n",
|
||||
"# NOTE: need this to force jupyter to reload imports:\n",
|
||||
"for key in list(sys.modules.keys()):\n",
|
||||
" if key.startswith('src.'):\n",
|
||||
" del sys.modules[key];\n",
|
||||
"\n",
|
||||
"os.chdir(os.path.dirname(_dh[0]));\n",
|
||||
"sys.path.insert(0, os.getcwd());\n",
|
||||
"\n",
|
||||
"from src.examples_dilations import *;"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# User input:\n",
|
||||
"N = 4; # dimension of the Hilbert space.\n",
|
||||
"d = 4;\n",
|
||||
"\n",
|
||||
"# If you ensure that the failure of S_{T,K} >= 0 only occurs for K = {1,2,...,d}\n",
|
||||
"# + you want S_TK > 0 (strictly) for all K ≠ {1,2,...,d}:\n",
|
||||
"alpha = 1/math.sqrt(d - 0.5);\n",
|
||||
"\n",
|
||||
"# If you ensure that the failure of S_{T,K} >= 0 only occurs for K = {1,2,...,d}\n",
|
||||
"# alpha = 1/math.sqrt(d - 1);\n",
|
||||
"\n",
|
||||
"# Otherwise:\n",
|
||||
"# alpha = 1;"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"\n",
|
||||
"The marginal semigroups T_i and their generators A_i:\n",
|
||||
"\n",
|
||||
" i spec bound of A_i A_i dissipative (<==> T_i contractive)?\n",
|
||||
"--- ------------------- -----------------------------------------\n",
|
||||
" 1 -0.465478 True\n",
|
||||
" 2 -0.465478 True\n",
|
||||
" 3 -0.465478 True\n",
|
||||
" 4 -0.465478 True\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# create the generators `A_i` of the marginal semigroups `T_i`:`\n",
|
||||
"A = [\n",
|
||||
" generate_semigroup_generator(\n",
|
||||
" shape = [N, N],\n",
|
||||
" rational = True,\n",
|
||||
" base = 100,\n",
|
||||
" alpha = alpha,\n",
|
||||
" )\n",
|
||||
" for _ in range(d)\n",
|
||||
"];\n",
|
||||
"\n",
|
||||
"data = [];\n",
|
||||
"for i, A_i in enumerate(A):\n",
|
||||
" omega_Re = spec_bounds((1/2)*(A_i + A_i.T.conj()));\n",
|
||||
" omega = spec_bounds(A_i);\n",
|
||||
" data.append((i+1, omega_Re, True if (omega_Re <= 0) else False));\n",
|
||||
"\n",
|
||||
"repr = tabulate(\n",
|
||||
" tabular_data = data,\n",
|
||||
" headers = ['i', 'spec bound of A_i', 'A_i dissipative (<==> T_i contractive)?'],\n",
|
||||
" showindex = False,\n",
|
||||
" floatfmt = '.6f',\n",
|
||||
" colalign = ['center', 'center', 'center'],\n",
|
||||
" tablefmt = 'simple',\n",
|
||||
");\n",
|
||||
"print(f'\\nThe marginal semigroups T_i and their generators A_i:\\n\\n{repr}');"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"\n",
|
||||
"Dissipation operators:\n",
|
||||
"\n",
|
||||
" K min σ(S_{T,K}) S_{T,K} >= 0?\n",
|
||||
"------------ ---------------- ---------------\n",
|
||||
" [] 1.000000 True\n",
|
||||
" [3] 0.465478 True\n",
|
||||
" [2] 0.465478 True\n",
|
||||
" [1] 0.465478 True\n",
|
||||
" [0] 0.465478 True\n",
|
||||
" [2, 3] 0.212675 True\n",
|
||||
" [1, 3] 0.183738 True\n",
|
||||
" [1, 2] 0.217453 True\n",
|
||||
" [0, 3] 0.200206 True\n",
|
||||
" [0, 2] 0.301332 True\n",
|
||||
" [0, 1] 0.215681 True\n",
|
||||
" [1, 2, 3] 0.058427 True\n",
|
||||
" [0, 2, 3] 0.075037 True\n",
|
||||
" [0, 1, 3] 0.056030 True\n",
|
||||
" [0, 1, 2] 0.077350 True\n",
|
||||
"[0, 1, 2, 3] -0.082403 False\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# compute the dissipation operators `S_TK` for each `K ⊆ {1,2,...,d}``:\n",
|
||||
"S, beta_T = dissipation_operators(shape=[N, N], A=A);\n",
|
||||
"\n",
|
||||
"data = [];\n",
|
||||
"for K, S_TK, b in sorted(S, key=lambda x: len(x[0])):\n",
|
||||
" data.append((K, b, True if b >= -MACHINE_EPS else False))\n",
|
||||
"\n",
|
||||
"repr = tabulate(\n",
|
||||
" tabular_data = data,\n",
|
||||
" headers = ['K', 'min σ(S_{T,K})', 'S_{T,K} >= 0?'],\n",
|
||||
" showindex = False,\n",
|
||||
" floatfmt = '.6f',\n",
|
||||
" colalign = ['center', 'center', 'center'],\n",
|
||||
" tablefmt = 'simple',\n",
|
||||
");\n",
|
||||
"print(f'\\nDissipation operators:\\n\\n{repr}');"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"\n",
|
||||
"β_T = min_K min σ(S_{T,K}) = -0.082403\n",
|
||||
"⟹ T is not compeletely dissipative.\n",
|
||||
"⟹ (by Thm 1.1) T does not have a regular unitary dilation.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Display summary:\n",
|
||||
"print('')\n",
|
||||
"print(f'β_T = min_K min σ(S_{{T,K}}) = {beta_T:.6f}');\n",
|
||||
"if beta_T == 0:\n",
|
||||
" print('⟹ T is compeletely dissipative.');\n",
|
||||
" print('⟹ (by Thm 1.1) T has a regular unitary dilation.');\n",
|
||||
"elif beta_T > 0:\n",
|
||||
" print('⟹ T is compeletely super dissipative.');\n",
|
||||
" print('⟹ (by Thm 1.1) T has a regular unitary dilation.');\n",
|
||||
"else:\n",
|
||||
" print('⟹ T is not compeletely dissipative.');\n",
|
||||
" print('⟹ (by Thm 1.1) T does not have a regular unitary dilation.');"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"jupytext": {
|
||||
"cell_metadata_filter": "-all"
|
||||
},
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3.10.6 64-bit",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.10.6"
|
||||
},
|
||||
"widgets": {
|
||||
"application/vnd.jupyter.widget-state+json": {
|
||||
"state": {},
|
||||
"version_major": 2,
|
||||
"version_minor": 0
|
||||
}
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 4
|
||||
}
|
||||
|
|
@ -0,0 +1,288 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
# AUTHOR: Raj Dahya
|
||||
# CREATED: 05.09.2022
|
||||
# EDITED: 08.10.2022
|
||||
# TYPE: Source code
|
||||
# DEPARTMENT: Fakult\"at for Mathematik und Informatik
|
||||
# INSTITUTE: Universit\"at Leipzig
|
||||
# NOTE: Collapse this block.
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
# IMPORTS
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
import re;
|
||||
import math
|
||||
from tkinter import E;
|
||||
import numpy as np;
|
||||
from typing import Generator;
|
||||
import itertools;
|
||||
from tabulate import tabulate;
|
||||
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
# CONSTANTS + SETTINGS
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
np.set_printoptions(precision=3);
|
||||
MACHINE_EPS = 0.5e-12;
|
||||
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
# SECONDARY METHODS - Generate C0-Semigroup Generator
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
def generate_semigroup_generator(
|
||||
shape: list[int, int],
|
||||
rational: bool = False,
|
||||
base: int = 2,
|
||||
alpha: float = 1,
|
||||
) -> np.ndarray:
|
||||
'''
|
||||
Creates bounded generators for d commuting C0-semigroups as in the paper (cf. §5.3).
|
||||
|
||||
@inputs
|
||||
- `shape` - <int, int> Desired shape of output.
|
||||
- `rational` - <bool> Whether or not entries are to be rational.
|
||||
- `base` - <int> If `rational = True`, fixes the denominator of the rational numbers.
|
||||
- `alpha` - <float> Additional parameter to scale the D_i operators.
|
||||
|
||||
NOTE: in the paper `α = 1` was chosen. However one can choose any value in `(1/√d, \infty)`.
|
||||
By choosing any value `α ≥ 1/√(d-1)`, by the computations in Proposition 5.3
|
||||
one can force that the S_TK operators only fail to be positive when |K| > d-1.
|
||||
|
||||
@returns
|
||||
Randomly created bounded generators for C0-semigroups as in the paper (cf. §5.3).
|
||||
'''
|
||||
sh = tuple(shape)
|
||||
m, n = sh;
|
||||
assert m == n, 'must be square';
|
||||
assert m % 2 == 0, 'dimensions must be divisible by 2';
|
||||
N = int(m/2);
|
||||
I = generate_op_identity(shape=[2*N, 2*N]);
|
||||
O = generate_op_zero(shape=[N, N]);
|
||||
U = generate_op_unitary(
|
||||
shape = [N, N],
|
||||
rational = rational,
|
||||
base = base,
|
||||
);
|
||||
D = alpha * U;
|
||||
A = -I + np.vstack([
|
||||
np.hstack([ O, -2*D ]),
|
||||
np.hstack([ O, O ]),
|
||||
]);
|
||||
return A;
|
||||
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
# SECONDARY METHODS - Compute spectral bounds
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
def spec_bounds(A: np.ndarray) -> float:
|
||||
'''
|
||||
Computes the spectral bounds of an operator:
|
||||
```tex
|
||||
\max_{i} \max \{ Re \lambda : \lambda \in \sigma(A_{i}) \}
|
||||
```
|
||||
'''
|
||||
sigma_A = np.linalg.eigvals(A);
|
||||
return max(np.real(sigma_A));
|
||||
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
# SECONDARY METHODS - Generate Dissipation Operator
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
def dissipation_operators(
|
||||
shape: list[int, int],
|
||||
A: list[np.ndarray],
|
||||
) -> tuple[
|
||||
list[tuple[
|
||||
list[int],
|
||||
np.ndarray,
|
||||
float,
|
||||
]],
|
||||
float,
|
||||
]:
|
||||
'''
|
||||
Computes the dissipation operators and their minimum spectral values.
|
||||
'''
|
||||
d = len(A);
|
||||
indexes = list(range(d));
|
||||
data = [];
|
||||
beta = np.inf;
|
||||
for K in power_set(indexes):
|
||||
S_TK, eig_K = dissipation_operator(shape=shape, A=A, K=K);
|
||||
data.append((K, S_TK, eig_K));
|
||||
beta = min(eig_K, beta);
|
||||
return data, beta;
|
||||
|
||||
def dissipation_operator(
|
||||
shape: list[int, int],
|
||||
A: list[np.ndarray],
|
||||
K: list[int],
|
||||
) -> tuple[np.ndarray, float]:
|
||||
'''
|
||||
Computes the dissipation operator `S_{T,K}` and its minimum spectral values.
|
||||
'''
|
||||
sh = tuple(shape)
|
||||
m, n = sh;
|
||||
assert m == n, 'must be square';
|
||||
N = m;
|
||||
O = generate_op_zero(shape=[N, N]);
|
||||
I = generate_op_identity(shape=[N, N]);
|
||||
S = O;
|
||||
for C, C_ in partitions_two(K):
|
||||
A_C = I;
|
||||
A_C_ = I;
|
||||
for i in C:
|
||||
A_C = A_C @ -A[i];
|
||||
for i in C_:
|
||||
A_C_ = A_C_ @ -A[i];
|
||||
S += A_C.conj().T @ A_C_;
|
||||
S = (1/2**len(K)) * S;
|
||||
eig = min(np.real(np.linalg.eigvals(S)));
|
||||
return S, eig;
|
||||
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
# MISCELLANEOUS METHODS - SETS
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
def power_set(K: list[int]) -> Generator[list[int], None, None]:
|
||||
'''
|
||||
Iterates through all subsets `C` of a set `K`.
|
||||
'''
|
||||
n = len(K);
|
||||
for ch in itertools.product(*([[0,1]]*n)):
|
||||
yield [ K[i] for i in range(n) if ch[i] == 1 ];
|
||||
|
||||
def partitions_two(K: list[int]) -> Generator[tuple[list[int],list[int]], None, None]:
|
||||
'''
|
||||
Iterates through all partitions `[C1, C2]` of a set `K` of size 2.
|
||||
'''
|
||||
n = len(K);
|
||||
for ch in itertools.product(*([[0,1]]*n)):
|
||||
yield (
|
||||
[ K[i] for i in range(n) if ch[i] == 1 ],
|
||||
[ K[i] for i in range(n) if ch[i] == 0 ],
|
||||
);
|
||||
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
# MISCELLANEOUS METHODS - BASIC OPERATORS
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
def generate_op_zero(shape: list[int, int]) -> np.ndarray:
|
||||
'''
|
||||
Returns the zero operator with the right shape and type.
|
||||
'''
|
||||
sh = tuple(shape)
|
||||
O = np.zeros(sh).astype(complex);
|
||||
return O;
|
||||
|
||||
def generate_op_identity(shape: list[int, int]) -> np.ndarray:
|
||||
'''
|
||||
Returns the identity operator with the right shape and type.
|
||||
'''
|
||||
sh = tuple(shape)
|
||||
m, n = sh;
|
||||
assert m == n, 'must be square';
|
||||
I = np.eye(m).astype(complex);
|
||||
return I;
|
||||
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
# MISCELLANEOUS METHODS - GENERATE RANDOM MATRICES
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
def generate_op_selfadjoint(
|
||||
shape: list[int, int],
|
||||
bound: float = 1.,
|
||||
) -> np.ndarray:
|
||||
'''
|
||||
Generates a self-adjoint operator.
|
||||
|
||||
@inputs
|
||||
- `shape` - <int, int> Desired shape of the output.
|
||||
- `bound` - <float> The desired norm of output operator.
|
||||
|
||||
@returns
|
||||
A randomly generated self-adjoint operator as a matrix.
|
||||
'''
|
||||
sh = tuple(shape)
|
||||
m, n = sh;
|
||||
assert m == n, 'must be square';
|
||||
A = np.random.random(sh) + 1j * np.random.random(sh);
|
||||
S = (A + A.conj().T)/2.;
|
||||
M = np.linalg.norm(S);
|
||||
if M > 0:
|
||||
S = (bound/M) * S;
|
||||
return S;
|
||||
|
||||
def generate_onb(shape: list[int, int]) -> np.ndarray:
|
||||
'''
|
||||
Generates an ONB.
|
||||
|
||||
@inputs
|
||||
- `shape` - <int, int> Desired shape of matrices.
|
||||
|
||||
@returns
|
||||
An ONB for an n-dim complex Hilbert space as a matrix.
|
||||
'''
|
||||
sh = tuple(shape)
|
||||
m, n = sh;
|
||||
assert m == n, 'must be square';
|
||||
N = m;
|
||||
O = generate_op_zero(shape=[N, N]);
|
||||
P = O;
|
||||
for k in range(N):
|
||||
while True:
|
||||
v = np.random.random((N,)) + 1j * np.random.random((N,));
|
||||
for i in range(k):
|
||||
# NOTE: np.inner does not work.
|
||||
coeff = np.sum([ v[j] * P[j, i].conj() for j in range(N)]);
|
||||
v -= coeff * P[:, i];
|
||||
norm_v = np.linalg.norm(v);
|
||||
if norm_v > 0.:
|
||||
P[:, k] = v / norm_v;
|
||||
break;
|
||||
return P;
|
||||
|
||||
def generate_op_unitary(
|
||||
shape: list[int, int],
|
||||
rational: bool = False,
|
||||
base: int = 2,
|
||||
) -> np.ndarray:
|
||||
'''
|
||||
Generates a unitary operator.
|
||||
|
||||
@inputs
|
||||
- `shape` - <int, int> Desired shape of output.
|
||||
- `rational` - <bool> Whether or not entries are to be rational.
|
||||
- `base` - <int> If `rational = True`, fixes the denominator of the rational numbers.
|
||||
|
||||
@returns
|
||||
A randomly generated unitary operator as a matrix.
|
||||
'''
|
||||
sh = tuple(shape)
|
||||
m, n = sh;
|
||||
assert m == n, 'must be square';
|
||||
N = m;
|
||||
P = generate_onb(shape=shape);
|
||||
if rational:
|
||||
theta = 2 * np.pi * np.random.randint(low=0, high=base, size=N) / base;
|
||||
else:
|
||||
theta = 2 * np.pi * np.random.random(size=N);
|
||||
D = np.diag(np.exp(1j * theta));
|
||||
U = P @ D @ P.conj().T;
|
||||
return U;
|
||||
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
# MISCELLANEOUS METHODS - Get user input
|
||||
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
def ask_confirmation(message: str):
|
||||
while True:
|
||||
answer = input(f'{message} (y/n): ');
|
||||
if re.match(pattern='^y(es)?|1$', string=answer, flags=re.IGNORECASE):
|
||||
return True;
|
||||
elif re.match(pattern='^n(o)?|0$', string=answer, flags=re.IGNORECASE):
|
||||
return False;
|
||||
Loading…
Reference in New Issue