Variational Quantum Singular Value Decomposition¶

import numpy as np
from numpy import pi as PI
from matplotlib import pyplot as plt
from scipy.stats import unitary_group
from scipy.linalg import norm

import paddle
from paddle_quantum.ansatz import Circuit
from paddle_quantum.linalg import dagger

# Generate matrix M
def M_generator():
    I = np.array([[1, 0], [0, 1]])
    Z = np.array([[1, 0], [0, -1]])
    X = np.array([[0, 1], [1, 0]])
    Y = np.array([[0, -1j], [1j, 0]])
    M = 2 *np.kron(X, Z) + 6 * np.kron(Z, X) + 3 * np.kron(I, I)
    return M.astype('complex64')

print('The matrix M we want to decompose is: ')
print(M_generator())

# We only need the following line of code to complete SVD
U, D, V_dagger = np.linalg.svd(M_generator(), full_matrices=True)


# Print decomposition results
print("The singular values of the matrix from large to small are:")
print(D)
print("The decomposed unitary matrix U is:")
print(U)
print("The decomposed unitary matrix V_dagger is:")
print(V_dagger)

# Then assemble it back, can we restore the original matrix?
M_reconst = np.matmul(U, np.matmul(np.diag(D), V_dagger))
print(M_reconst)

# First fix the random seed, in order to reproduce the results at any time
np.random.seed(42)

# Set the number of qubits, which determines the dimension of the Hilbert space
N = 3

# Make a random matrix generator
def random_M_generator():
    M = np.random.randint(10, size = (2**N, 2**N)) + 1j*np.random.randint(10, size = (2**N, 2**N))
    return M

M = random_M_generator()
M_err = np.copy(M)


# Output the matrix M
print('The matrix M we want to decompose is:')
print(M)

# Apply SVD and record the exact singular values
U, D, V_dagger = np.linalg.svd(M, full_matrices=True)
print("The singular values of the matrix M are:")
print(D)

# Hyperparameter settings
RANK = 8       # Set the number of rank you want to learn
ITR = 100   # Number of iterations
LR = 0.02   # Learning rate
SEED = 14   # Random seed

# Set the learning weight 
weight = np.arange(3 * RANK, 0, -3).astype('complex128')
print('The selected weight is:')
print(weight)

# Record the optimization process
loss_list, singular_value_list = [], []
U_learned, V_dagger_learned = [], []

# Construct the VQSVD network and train
net = VQSVD(matrix=M, weights=weight, num_qubits=num_qubits, depth=cir_depth, rank=RANK, lr=LR, itr=ITR, seed=SEED)
loss_list, singular_value_list = net.train()

# Draw a learning curve
loss_plot(loss_list)

# Record the last two learned unitary matrices
U_learned = net.get_matrix_U().numpy()
V_dagger_learned = dagger(net.get_matrix_V()).numpy()

singular_value = singular_value_list[-1]
err_subfull, err_local, err_SVD = [], [], []
U, D, V_dagger = np.linalg.svd(M, full_matrices=True)


# Calculate the Frobenius-norm error
for i in range(RANK):
    lowrank_mat = np.matrix(U[:, :i]) * np.diag(D[:i])* np.matrix(V_dagger[:i, :])
    recons_mat = np.matrix(U_learned[:, :i]) * np.diag(singular_value[:i])* np.matrix(V_dagger_learned[:i, :])
    err_local.append(norm(lowrank_mat - recons_mat)) 
    err_subfull.append(norm(M_err - recons_mat))
    err_SVD.append(norm(M_err- lowrank_mat))


# Plot
fig, ax = plt.subplots()
ax.plot(list(range(1, RANK+1)), err_subfull, "o-.", 
        label = 'Reconstruction via VQSVD')
ax.plot(list(range(1, RANK+1)), err_SVD, "^--", 
        label='Reconstruction via SVD')
plt.xlabel('Singular Value Used (Rank)', fontsize = 14)
plt.ylabel('Norm Distance', fontsize = 14)
leg = plt.legend(frameon=True)
leg.get_frame().set_edgecolor('k')

# Then we look at the effect of the classic singular value decomposition
U, sigma, V = np.linalg.svd(imgmat)

for i in range(5, 16, 5):
    reconstimg = np.matrix(U[:, :i]) * np.diag(sigma[:i]) * np.matrix(V[:i, :])
    plt.imshow(reconstimg, cmap='gray')
    title = f"n = {i}"
    plt.title(title)
    plt.show()

# Record the optimization process
loss_list, singular_value_list = [], []
U_learned, V_dagger_learned = [], []
    
# Construct the VQSVD network and train
net = VQSVD(matrix=Mat_generator(), weights=weight, num_qubits=num_qubits, depth=cir_depth, rank=RANK, lr=LR, itr=ITR, seed=SEED)
loss_list, singular_value_list = net.train()

# Record the last two unitary matrices learned
U_learned = net.get_matrix_U().numpy()
V_dagger_learned = dagger(net.get_matrix_V()).numpy()

singular_value = singular_value_list[-1]
mat = np.matrix(U_learned.real[:, :RANK]) * np.diag(singular_value[:RANK])* np.matrix(V_dagger_learned.real[:RANK, :])

reconstimg = mat
plt.imshow(reconstimg, cmap='gray')