Graphical illustration

Code

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import wishart
import pandas as pd
import plotly.graph_objects as go
import time
import pickle
from my_functions import Ellipses
from my_functions import Distance
from my_functions import mat_to_vec
from my_functions import ClusterMatricesCPC
from my_functions import ClusterMatrices
from my_functions import ClusteringCPC
from sklearn.cluster import KMeans
from ipywidgets import IntProgress
from IPython.display import display
from time import sleep
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.io as pio
from tqdm import tqdm

Code

# functions
def rotation_matrix(theta):
    return np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])


def equal_diag(sigma):
    thetas = np.linspace(0, np.pi, 1000)
    rot_mats = np.transpose(rotation_matrix(thetas), (2, 0, 1))
    rotated_Psis = rot_mats @ sigma @ np.transpose(rot_mats, (0, 2, 1))
    diff = np.array(
        list(map(lambda x: np.square((np.diag(x)[0] - np.diag(x)[1])), rotated_Psis))
    )
    indx = diff.argmin()
    theta = thetas[indx]
    rot_mat = rotation_matrix(theta)
    sigma_tran = rot_mat @ sigma @ rot_mat.T
    return rot_mat, sigma_tran


def cpca_sigma(sigma, Psis):
    _, eigenvectors = np.linalg.eig(sigma)
    eigenvalues = (eigenvectors.T @ Psis @ eigenvectors) * np.eye(2)
    Psishat = eigenvectors @ eigenvalues @ eigenvectors.T
    return Psishat


def principal_tensors(sigma, Psis):
    _, eigenvectors = np.linalg.eig(sigma)
    eigenvalues = (eigenvectors.T @ Psis @ eigenvectors) * np.eye(2)
    Psishat = eigenvectors @ eigenvalues @ eigenvectors.T
    Xshat = mat_to_vec(Psishat)
    return Psishat, Xshat, eigenvectors, eigenvalues


# change a vector to a upper triangular matrix
def vec_to_mat(Xs):
    mat = np.diag(Xs[[0, 2]])
    mat[0, 1] = mat[1, 0] = Xs[1]
    return mat

Simulation 1

Code

# simulation
theta1 = np.pi / 6
theta2 = theta1 + np.pi / 15
df = 20
n = 50000

eigenvectors1 = np.array(
    [[np.cos(theta1), -np.sin(theta1)], [np.sin(theta1), np.cos(theta1)]]
)
eigenvectors2 = np.array(
    [[np.cos(theta2), -np.sin(theta2)], [np.sin(theta2), np.cos(theta2)]]
)

# eigenvalues1 = np.diag([2, 0.5])
# eigenvalues2 = np.diag([3, 0.9])

# sigma1 = eigenvectors1 @ eigenvalues1 @ eigenvectors1.T
sigma1 = np.array([[1, 0.75], [0.75, 1]])
# sigma1 = np.array([[1, 0.999999], [0.999999, 1]])
Psis1 = wishart.rvs(df=df, scale=sigma1, size=int(n)) / df
Psis2 = np.array(list(map(lambda x: np.cov(x), Psis1)))
# Psis1 = np.array(list(map(lambda x: np.linalg.inv(x), Psis1)))


Xs1 = mat_to_vec(Psis1)
Xs2 = mat_to_vec(Psis2)

\[\begin{align} \mathbf \Psi \sim \mathcal W_2 \left(\begin{bmatrix} 1 & 0.25 \\ 0.25 & 1 \end{bmatrix}, 40\right ) \end{align}\] where the diagonal of the covariance matrix are the same value.

Visulize the data \(\mathop{\mathrm{vec}}{\left(\mathbf \Psi\right)}\)

Code

fig = go.Figure(
    data=[
        go.Scatter3d(
            x=Xs1[:, 0],
            y=Xs1[:, 1],
            z=Xs1[:, 2],
            mode="markers",
            marker=dict(size=2, color="gray", opacity=0.5),
        )
    ]
)
# add Xs2
# fig.add_trace(
#     go.Scatter3d(
#         x=Xs2[:, 0],
#         y=Xs2[:, 1],
#         z=Xs2[:, 2],
#         mode="markers",
#         marker=dict(size=2, color="red"),
#     )
# )
fig.update_layout(
    width=800,
    height=800,
    scene=dict(
        xaxis=dict(title="X1"),
        yaxis=dict(title="X2"),
        zaxis=dict(title="X3"),
    ),
)
# fig.update_traces(marker=dict(size=2, color="gray", opacity=0.5))
fig.update_layout(
    scene=dict(
        # aspectmode="cube",
        xaxis=dict(backgroundcolor="white", gridcolor="gray"),
        yaxis=dict(backgroundcolor="white", gridcolor="gray"),
        zaxis=dict(backgroundcolor="white", gridcolor="gray"),
    )
)

Visulize the principal positive semi-definite tensors

\[\begin{align} \mathcal P_{\mathbf B}(\mathbf \Psi) = \mathbf B\left(\mathbf B^T \mathbf \Psi \mathbf B \circ \mathbf I\right) \mathbf B^T \end{align}\]

Code

# plot Xs in 3d
# fig size
fig = go.Figure(
    data=[
        go.Scatter3d(
            x=Xs1[:, 0],
            y=Xs1[:, 1],
            z=Xs1[:, 2],
            mode="markers",
            marker=dict(size=2, color="blue"),
        )
    ]
)
fig.update_layout(
    width=800,
    height=800,
    scene=dict(
        xaxis=dict(title="X1"),
        yaxis=dict(title="X2"),
        zaxis=dict(title="X3"),
    ),
)
fig.update_traces(marker=dict(size=2, color="gray", opacity=0.5))

Psishat, Xshat, eigenvectors, eigenvalues = principal_tensors(sigma1, Psis1)
# add Xshat as a surface

# sort Xshat by the distance to the origin for three columns
indx = np.square((Xshat - Xshat.mean(axis=0))).sum(axis=1).argmin()
Xshat = np.concatenate(
    (Xshat[indx:, :], Xshat[:indx, :]), axis=0
)  # move the point closest to the origin to the first row

# sqrtn = int(np.sqrt(n))
# fig.add_trace(
#     go.Surface(
#         x=Xshat[:, 0].reshape(sqrtn, sqrtn),
#         y=Xshat[:, 1].reshape(sqrtn, sqrtn),
#         z=Xshat[:, 2].reshape(sqrtn, sqrtn),
#         opacity=0.2,
#         colorscale="Viridis",
#     )
# )


fig.add_trace(
    go.Scatter3d(
        x=Xshat[:, 0],
        y=Xshat[:, 1],
        z=Xshat[:, 2],
        mode="markers",
        marker=dict(size=2, color="black", opacity=0.2),
    )
)
indx = np.square((Xshat - Xshat.mean(axis=0))).sum(axis=1).argmin()

# add the point Xshat[indx] on figure
# fig.add_trace(
#     go.Scatter3d(
#         x=[Xshat[indx, 0]],
#         y=[Xshat[indx, 1]],
#         z=[Xshat[indx, 2]],
#         mode="markers",
#         marker=dict(size=8, color="red"),
#     )
# )


def off_diag(x):
    mat = np.eye(2) - np.eye(2)
    mat[0, 1] = mat[1, 0] = x
    return mat


eigenvalue_line = (
    np.array(list(map(lambda x: off_diag(x), np.linspace(-1, 1, 1000))))
    + eigenvalues[indx]
)

Psis_line = eigenvectors @ eigenvalue_line @ eigenvectors.T

Xline = mat_to_vec(Psis_line)


# add the line on figure
fig.add_trace(
    go.Scatter3d(
        x=Xline[:, 0],
        y=Xline[:, 1],
        z=Xline[:, 2],
        mode="lines",
        line=dict(width=10, color="red"),
    )
)

fig.update_layout(
    showlegend=False,
    scene=dict(
        xaxis=dict(title="X1"),
        yaxis=dict(title="X2"),
        zaxis=dict(title="X3"),
    ),
    width=800,
    height=800,
)

# background color to white and grid color to gray

fig.update_layout(
    scene=dict(
        # aspectmode="cube",
        xaxis=dict(backgroundcolor="white", gridcolor="gray"),
        yaxis=dict(backgroundcolor="white", gridcolor="gray"),
        zaxis=dict(backgroundcolor="white", gridcolor="gray"),
    )
)
# fig.update_layout(
#     scene=dict(
#         xaxis=dict(backgroundcolor="white", gridcolor="white"),
#         yaxis=dict(backgroundcolor="white", gridcolor="white"),
#         zaxis=dict(backgroundcolor="white", gridcolor="white"),
#     )
# )

# remove x y z axis
# fig.update_layout(
#     scene=dict(
#         xaxis=dict(showticklabels=False, visible=False),
#         yaxis=dict(showticklabels=False, visible=False),
#         zaxis=dict(showticklabels=False, visible=False),
#     )
# )

fig.show()

Simulation 2

unequal diagonal elements

Code

aplha = 4
sigma1 = np.array([[aplha, 1], [1, 1.5]])
Psis1 = wishart.rvs(df=df, scale=sigma1, size=int(n)) / df
Xs1 = mat_to_vec(Psis1)

\[\begin{align} \mathbf \Psi \sim \mathcal W_2 \left(\begin{bmatrix} 4 & 1 \\ 1 & 1.5 \end{bmatrix}, 40\right ) \end{align}\] where the diagonal of the covariance matrix are the same value.

Visulize the data \(\mathop{\mathrm{vec}}{\left(\mathbf \Psi\right)}\)

Code

fig = go.Figure(
    data=[
        go.Scatter3d(
            x=Xs1[:, 0],
            y=Xs1[:, 1],
            z=Xs1[:, 2],
            mode="markers",
            marker=dict(size=2, color="blue"),
        )
    ]
)
fig.update_layout(
    width=800,
    height=800,
    scene=dict(
        xaxis=dict(title="X1"),
        yaxis=dict(title="X2"),
        zaxis=dict(title="X3"),
    ),
)
fig.update_traces(marker=dict(size=2, color="gray", opacity=0.5))
fig.update_layout(
    scene=dict(
        aspectmode="cube",
        xaxis=dict(backgroundcolor="white", gridcolor="gray", range=[0, 2 * aplha]),
        yaxis=dict(backgroundcolor="white", gridcolor="gray", range=[-aplha, aplha]),
        zaxis=dict(backgroundcolor="white", gridcolor="gray", range=[-aplha, aplha]),
    )
)

Visulize the principal positive semi-definite tensors

\[\begin{align} \mathcal P_{\mathbf B}(\mathbf \Psi) = \mathbf B\left(\mathbf B^T \mathbf \Psi \mathbf B \circ \mathbf I\right) \mathbf B^T \end{align}\]

Code

# plot Xs in 3d
# fig size
fig = go.Figure(
    data=[
        go.Scatter3d(
            x=Xs1[:, 0],
            y=Xs1[:, 1],
            z=Xs1[:, 2],
            mode="markers",
            marker=dict(size=2, color="blue"),
        )
    ]
)
fig.update_layout(
    width=800,
    height=800,
    scene=dict(
        xaxis=dict(title="X1"),
        yaxis=dict(title="X2"),
        zaxis=dict(title="X3"),
    ),
)
fig.update_traces(marker=dict(size=2, color="gray", opacity=0.5))

Psishat, Xshat, eigenvectors, eigenvalues = principal_tensors(sigma1, Psis1)
# add Xshat as a surface

# sort Xshat by the distance to the origin for three columns
indx = np.square((Xshat - Xshat.mean(axis=0))).sum(axis=1).argmin()
Xshat = np.concatenate(
    (Xshat[indx:, :], Xshat[:indx, :]), axis=0
)  # move the point closest to the origin to the first row

# sqrtn = int(np.sqrt(n))
# fig.add_trace(
#     go.Surface(
#         x=Xshat[:, 0].reshape(sqrtn, sqrtn),
#         y=Xshat[:, 1].reshape(sqrtn, sqrtn),
#         z=Xshat[:, 2].reshape(sqrtn, sqrtn),
#         opacity=0.2,
#         colorscale="Viridis",
#     )
# )

fig.add_trace(
    go.Scatter3d(
        x=Xshat[:, 0],
        y=Xshat[:, 1],
        z=Xshat[:, 2],
        mode="markers",
        marker=dict(size=2, color="black", opacity=0.2),
    )
)
indx = np.square((Xshat - Xshat.mean(axis=0))).sum(axis=1).argmin()

# add the point Xshat[indx] on figure
# fig.add_trace(
#     go.Scatter3d(
#         x=[Xshat[indx, 0]],
#         y=[Xshat[indx, 1]],
#         z=[Xshat[indx, 2]],
#         mode="markers",
#         marker=dict(size=8, color="red"),
#     )
# )


def off_diag(x):
    mat = np.eye(2) - np.eye(2)
    mat[0, 1] = mat[1, 0] = x
    return mat


eigenvalue_line = (
    np.array(list(map(lambda x: off_diag(x), np.linspace(-1, 1, 1000))))
    + eigenvalues[indx]
)

Psis_line = eigenvectors @ eigenvalue_line @ eigenvectors.T

Xline = mat_to_vec(Psis_line)


# add the line on figure
fig.add_trace(
    go.Scatter3d(
        x=Xline[:, 0],
        y=Xline[:, 1],
        z=Xline[:, 2],
        mode="lines",
        line=dict(width=10, color="red"),
    )
)


# background color to white and grid color to gray

fig.update_layout(
    scene=dict(
        aspectratio=dict(x=1, y=1, z=1),
        xaxis=dict(backgroundcolor="white", gridcolor="gray", range=[0, 2 * aplha]),
        yaxis=dict(backgroundcolor="white", gridcolor="gray", range=[-aplha, aplha]),
        zaxis=dict(backgroundcolor="white", gridcolor="gray", range=[-aplha, aplha]),
    )
)

fig.update_layout(
    showlegend=False,
    scene=dict(
        xaxis=dict(title="X1"),
        yaxis=dict(title="X2"),
        zaxis=dict(title="X3"),
    ),
    width=800,
    height=800,
)
# fig.update_layout(
#     scene=dict(
#         xaxis=dict(backgroundcolor="white", gridcolor="white"),
#         yaxis=dict(backgroundcolor="white", gridcolor="white"),
#         zaxis=dict(backgroundcolor="white", gridcolor="white"),
#     )
# )

# remove x y z axis
# fig.update_layout(
#     scene=dict(
#         xaxis=dict(showticklabels=False, visible=False),
#         yaxis=dict(showticklabels=False, visible=False),
#         zaxis=dict(showticklabels=False, visible=False),
#     )
# )

fig.show()

Decomposition of \(\mathbf \Psi\)

find a orthonormal matrix \(\mathbf A\) s.t. \[\begin{align} \left( \mathbf A^T \mathbf \Psi \mathbf A \right) \circ \mathbf I - \frac{1}{p} \left[\mathbf 1^T \left(\left( \mathbf A^T \mathbf \Psi \mathbf A \right) \circ \mathbf I \right)\mathbf 1 \right] \mathbf I = \mathbf 0 \end{align}\]

Code

eigenvalues_Psis1 = np.array(list(map(lambda x: np.max(np.linalg.eig(x)[0]), Psis1)))

Code

eigenvalues_Psishat = np.array(list(map(lambda x: np.max(np.linalg.eig(x)[0]), Psishat)))

Code

# use density instead of frequency

plt.figure(figsize=(8, 5))
plt.hist(eigenvalues_Psis1, bins = 50, alpha=0.2, density = True, color = 'red')
plt.hist(eigenvalues_Psishat, bins = 50, alpha=0.2, density = True, color = 'blue')
plt.show()