MatrixAnalysis.hpp

Go to the documentation of this file.

// Copyright Quantinuum
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
 
#pragma once
 
#include <array>
#include <tklog/TketLog.hpp>
#include <vector>
 
#include "EigenConfig.hpp"
#include "UnitID.hpp"
#include "tket/Utils/Constants.hpp"
#include "tket/Utils/Expression.hpp"
 
namespace tket {
 
typedef Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic> MatrixXb;
typedef Eigen::Matrix<bool, Eigen::Dynamic, 1> VectorXb;
typedef Eigen::Matrix<bool, 2, 1> Vector2b;
// Eigen::Matrix8cd doesn't exist.
typedef Eigen::Matrix<std::complex<double>, 8, 8> Matrix8cd;
 
bool is_unitary(const Eigen::MatrixXcd &U, double tol = EPS);
bool is_projector(const Eigen::MatrixXcd &P, double tol = EPS);
 
Eigen::PermutationMatrix<Eigen::Dynamic> lift_perm(
    const std::map<unsigned, unsigned> &p);
 
// Convert between ILO-BE and DLO-BE conventions
Eigen::Matrix4cd reverse_indexing(const Eigen::Matrix4cd &m);
Matrix8cd reverse_indexing(const Matrix8cd &m);
Eigen::MatrixXcd reverse_indexing(const Eigen::MatrixXcd &m);
Eigen::VectorXcd reverse_indexing(const Eigen::VectorXcd &v);
 
Eigen::MatrixXcd apply_qubit_permutation(
    const Eigen::MatrixXcd &m, const qubit_map_t &perm);
Eigen::VectorXcd apply_qubit_permutation(
    const Eigen::VectorXcd &v, const qubit_map_t &perm);
 
std::pair<MatrixXb, MatrixXb> binary_LLT_decomposition(const MatrixXb &a);
std::vector<std::pair<unsigned, unsigned>> gaussian_elimination_col_ops(
    const MatrixXb &a, unsigned blocksize = 6);
std::vector<std::pair<unsigned, unsigned>> gaussian_elimination_row_ops(
    const MatrixXb &a, unsigned blocksize = 6);
 
std::tuple<Eigen::Matrix4cd, std::array<double, 3>, Eigen::Matrix4cd>
get_information_content(const Eigen::Matrix4cd &X);
 
// given a 4x4 unitary matrix (ILO-BE), returns two 2x2 unitaries that
// approximately make the input by kronecker product
std::pair<Eigen::Matrix2cd, Eigen::Matrix2cd> kronecker_decomposition(
    Eigen::Matrix4cd &U);
 
unsigned get_matrix_size(unsigned number_of_qubits);
 
unsigned get_number_of_qubits(unsigned matrix_size);
 
typedef Eigen::Triplet<std::complex<double>> TripletCd;
 
typedef Eigen::SparseMatrix<std::complex<double>> SparseMatrixXcd;
 
SparseMatrixXcd get_sparse_matrix(
    const std::vector<TripletCd> &triplets, unsigned rows, unsigned cols);
 
SparseMatrixXcd get_sparse_square_matrix(
    const std::vector<TripletCd> &triplets, unsigned rows);
 
std::vector<TripletCd> get_triplets(
    const SparseMatrixXcd &matr, double abs_epsilon = EPS);
 
std::vector<TripletCd> get_triplets(
    const Eigen::MatrixXcd &matr, double abs_epsilon = EPS);
 
double trace_fidelity(double a, double b, double c);
 
bool in_weyl_chamber(const std::array<Expr, 3> &k);
 
Eigen::Matrix2cd nth_root(const Eigen::Matrix2cd &u, unsigned long long n);
 
template <class MatrixT>
static inline MatrixT clamp_to_unitary(const MatrixT &A) {
  if (is_unitary(A)) {
    return A;
  } else {
    tket_log()->warn(
        "Non-unitary product of matrices assumed unitary: "
        "presuming rounding error and applying correction.");
    // Fact: if A is an arbitrary complex matrix, and
    //   A = U * S * V^+
    // is its singular-value decomposition, then
    //   U * V^+
    // is the unitary matrix closest to A in the Frobenius norm.
    // See e.g. https://math.stackexchange.com/a/2215371 for a proof.
    Eigen::JacobiSVD<MatrixT, Eigen::NoQRPreconditioner> svd(
        A, Eigen::DecompositionOptions::ComputeFullU |
               Eigen::DecompositionOptions::ComputeFullV);
    return svd.matrixU() * svd.matrixV().adjoint();
  }
}
 
template <class MatrixT>
MatrixT unitary_product2(const MatrixT &U, const MatrixT &V) {
  MatrixT A = U * V;
  return clamp_to_unitary(A);
}
 
template <class MatrixT>
MatrixT unitary_product3(const MatrixT &U, const MatrixT &V, const MatrixT &W) {
  MatrixT A = U * V * W;
  return clamp_to_unitary(A);
}
 
}  // namespace tket