✨ Enhance Matrix functionality with new methods and optimizations (#1802)

Author: simon1hofmann

CHANGELOG.md

@@ -42,7 +42,7 @@ with the exception that minor releases may include breaking changes.
   [#1569], [#1570], [#1572], [#1573], [#1580], [#1602], [#1620], [#1623],
   [#1624], [#1626], [#1627], [#1635], [#1638], [#1673], [#1675], [#1700],
   [#1710], [#1717], [#1728], [#1730], [#1749], [#1751], [#1762], [#1765],
-  [#1774], [#1780], [#1781], [#1782], [#1787])
+  [#1774], [#1780], [#1781], [#1782], [#1787], [#1802])
   ([**@burgholzer**], [**@denialhaag**], [**@taminob**], [**@DRovara**],
   [**@li-mingbao**], [**@Ectras**], [**@MatthiasReumann**],
   [**@simon1hofmann**])
@@ -598,6 +598,7 @@ changelogs._
 
 <!-- PR links -->
 
+[#1802]: https://github.com/munich-quantum-toolkit/core/pull/1802
 [#1787]: https://github.com/munich-quantum-toolkit/core/pull/1787
 [#1782]: https://github.com/munich-quantum-toolkit/core/pull/1782
 [#1781]: https://github.com/munich-quantum-toolkit/core/pull/1781

mlir/include/mlir/Dialect/QCO/Utils/Matrix.h

@@ -10,6 +10,9 @@
 
 #pragma once
 
+#include <llvm/ADT/ArrayRef.h>
+#include <mlir/Support/LLVM.h>
+
 #include <array>
 #include <complex>
 #include <cstddef>
@@ -26,6 +29,8 @@ using Complex = std::complex<double>;
 inline constexpr double MATRIX_TOLERANCE = 1e-14;
 
 class DynamicMatrix;
+struct Matrix4x4;
+struct SymmetricEigen4;
 
 /**
  * @brief 1x1 matrix for global-phase gates.
@@ -183,6 +188,12 @@ struct Matrix2x2 {
    */
   [[nodiscard]] Matrix2x2 adjoint() const;
 
+  /**
+   * @brief Returns the (non-conjugate) transpose of this matrix.
+   * @return Transposed matrix `A^T`.
+   */
+  [[nodiscard]] Matrix2x2 transpose() const;
+
   /**
    * @brief Returns the trace of this matrix.
    * @return Sum of diagonal entries.
@@ -195,6 +206,13 @@ struct Matrix2x2 {
    */
   [[nodiscard]] Complex determinant() const;
 
+  /**
+   * @brief Checks whether this matrix is approximately the identity.
+   * @param tol Maximum allowed complex modulus of each entry difference.
+   * @return True if every entry is within @p tol of the identity.
+   */
+  [[nodiscard]] bool isIdentity(double tol = MATRIX_TOLERANCE) const;
+
   /**
    * @brief Checks approximate equality using an absolute entry-wise tolerance.
    *
@@ -215,6 +233,30 @@ struct Matrix2x2 {
    * @return `true` when @p src is 2x2.
    */
   [[nodiscard]] bool assignFrom(const DynamicMatrix& src);
+
+  /**
+   * @brief Embed this single-qubit matrix into an @p numQubits-qubit Hilbert
+   * space.
+   *
+   * Wire @p qubitIndex uses the same MSB-first convention as @ref
+   * Matrix4x4::kron (high bit first operand, low bit second). For each basis
+   * pair whose untouched wires match, copies this matrix at the target qubit's
+   * row/column bits.
+   *
+   * @param numQubits Number of qubits in the target Hilbert space.
+   * @param qubitIndex Wire index to act on.
+   * @return Embedded unitary as a dynamic matrix.
+   */
+  [[nodiscard]] DynamicMatrix embedInNqubit(std::size_t numQubits,
+                                            std::size_t qubitIndex) const;
+
+  /**
+   * @brief Embed this single-qubit matrix into a two-qubit Hilbert space.
+   *
+   * @param qubitIndex Wire index (`0` = high bit / MSB, `1` = low bit).
+   * @return The `4x4` embedded unitary.
+   */
+  [[nodiscard]] Matrix4x4 embedInTwoQubit(std::size_t qubitIndex) const;
 };
 
 /**
@@ -322,6 +364,12 @@ struct Matrix4x4 {
    */
   [[nodiscard]] Matrix4x4 adjoint() const;
 
+  /**
+   * @brief Returns the (non-conjugate) transpose of this matrix.
+   * @return Transposed matrix `A^T`.
+   */
+  [[nodiscard]] Matrix4x4 transpose() const;
+
   /**
    * @brief Returns the trace of this matrix.
    * @return Sum of diagonal entries.
@@ -334,6 +382,81 @@ struct Matrix4x4 {
    */
   [[nodiscard]] Complex determinant() const;
 
+  /**
+   * @brief Checks whether this matrix is approximately the identity.
+   * @param tol Maximum allowed complex modulus of each entry difference.
+   * @return True if every entry is within @p tol of the identity.
+   */
+  [[nodiscard]] bool isIdentity(double tol = MATRIX_TOLERANCE) const;
+
+  /**
+   * @brief Returns the four diagonal entries `(m00, m11, m22, m33)`.
+   * @return Array of diagonal entries.
+   */
+  [[nodiscard]] std::array<Complex, K_ROWS> diagonal() const;
+
+  /**
+   * @brief Builds a diagonal matrix from four diagonal entries.
+   * @param diagonalEntries Diagonal entries `(m00, m11, m22, m33)`; must have
+   *        length `K_ROWS`.
+   * @return Diagonal matrix with the given entries.
+   */
+  [[nodiscard]] static Matrix4x4
+  fromDiagonal(ArrayRef<Complex> diagonalEntries);
+
+  /**
+   * @brief Kronecker product `lhs (x) rhs` of two single-qubit matrices.
+   *
+   * Uses the computational-basis bit order where the first operand labels the
+   * high bit, matching `UnitaryOpInterface::getUnitaryMatrix4x4`.
+   *
+   * @param lhs Left factor (acts on the high bit / qubit 0).
+   * @param rhs Right factor (acts on the low bit / qubit 1).
+   * @return The `4x4` Kronecker product.
+   */
+  [[nodiscard]] static Matrix4x4 kron(const Matrix2x2& lhs,
+                                      const Matrix2x2& rhs);
+
+  /**
+   * @brief Returns the entries of column @p col, top to bottom.
+   * @param col Column index in `[0, K_COLS)`.
+   * @return Array of the four column entries.
+   */
+  [[nodiscard]] std::array<Complex, K_ROWS> column(std::size_t col) const;
+
+  /**
+   * @brief Overwrites column @p col with @p values.
+   * @param col Column index in `[0, K_COLS)`.
+   * @param values New column entries, top to bottom; must have length `K_ROWS`.
+   */
+  void setColumn(std::size_t col, ArrayRef<Complex> values);
+
+  /**
+   * @brief Returns the entries of row @p row, left to right.
+   * @param row Row index in `[0, K_ROWS)`.
+   * @return View over the four row entries.
+   */
+  [[nodiscard]] ArrayRef<const Complex> row(std::size_t row) const;
+
+  /**
+   * @brief Overwrites row @p row with @p values.
+   * @param row Row index in `[0, K_ROWS)`.
+   * @param values New row entries, left to right; must have length `K_COLS`.
+   */
+  void setRow(std::size_t row, ArrayRef<Complex> values);
+
+  /**
+   * @brief Returns the element-wise real parts in row-major order.
+   * @return Real parts of all entries.
+   */
+  [[nodiscard]] std::array<double, K_SIZE_AT_COMPILE_TIME> realPart() const;
+
+  /**
+   * @brief Returns the element-wise imaginary parts in row-major order.
+   * @return Imaginary parts of all entries.
+   */
+  [[nodiscard]] std::array<double, K_SIZE_AT_COMPILE_TIME> imagPart() const;
+
   /**
    * @brief Checks approximate equality using an absolute entry-wise tolerance.
    *
@@ -354,6 +477,58 @@ struct Matrix4x4 {
    * @return `true` when @p src is 4x4.
    */
   [[nodiscard]] bool assignFrom(const DynamicMatrix& src);
+
+  /**
+   * @brief Embed this two-qubit matrix into an @p numQubits-qubit Hilbert
+   * space.
+   *
+   * Operand 0 labels the high bit of the pair and acts on @p q0Index; operand 1
+   * labels the low bit and acts on @p q1Index. For each basis pair whose other
+   * wires match, copies this matrix at the packed two-qubit row/column indices.
+   *
+   * @param numQubits Number of qubits in the target Hilbert space.
+   * @param q0Index Wire index of operand 0.
+   * @param q1Index Wire index of operand 1.
+   * @return Embedded unitary as a dynamic matrix.
+   */
+  [[nodiscard]] DynamicMatrix embedInNqubit(std::size_t numQubits,
+                                            std::size_t q0Index,
+                                            std::size_t q1Index) const;
+
+  /**
+   * @brief Reorder this matrix to act on qubits `{0, 1}`.
+   *
+   * @param q0Index Wire index of operand 0; @p q1Index wire index of operand 1.
+   * @return Reordered copy of this matrix.
+   */
+  [[nodiscard]] Matrix4x4 reorderForQubits(std::size_t q0Index,
+                                           std::size_t q1Index) const;
+
+  /**
+   * @brief Computes the eigendecomposition of this real symmetric matrix.
+   *
+   * @copydoc Matrix4x4::symmetricEigen4(const std::array<double, 16>&)
+   *
+   * @pre Entries are real (imaginary parts must be negligible). The real parts
+   * must form a symmetric matrix; imaginary parts are ignored.
+   */
+  [[nodiscard]] SymmetricEigen4 symmetricEigen4() const;
+
+  /**
+   * @brief Computes the eigendecomposition of a real symmetric `4x4` matrix.
+   *
+   * Uses Householder tridiagonalization (EISPACK `tred2`) followed by implicit
+   * QL iteration (`tql2`) on the tridiagonal form.
+   *
+   * @pre @p symmetric is real and symmetric: `symmetric[i,j] == symmetric[j,i]`
+   * for all `i, j`. Only the lower triangle (including the diagonal) is read,
+   * but supplying a non-symmetric matrix yields undefined numerical results.
+   *
+   * @param symmetric Row-major real symmetric `4x4` matrix.
+   * @return Ascending eigenvalues and matching eigenvectors (as columns).
+   */
+  [[nodiscard]] static SymmetricEigen4
+  symmetricEigen4(const std::array<double, 16>& symmetric);
 };
 
 /**
@@ -539,6 +714,41 @@ class DynamicMatrix {
   [[nodiscard]] bool isApprox(const DynamicMatrix& other,
                               double tol = MATRIX_TOLERANCE) const;
 
+  /**
+   * @brief Returns the trace of this matrix.
+   * @return Sum of diagonal entries.
+   */
+  [[nodiscard]] Complex trace() const;
+
+  /**
+   * @brief Matrix product `*this * rhs`.
+   * @param rhs Right-hand factor.
+   * @return Product of the two matrices.
+   */
+  [[nodiscard]] DynamicMatrix operator*(const DynamicMatrix& rhs) const;
+
+  /**
+   * @brief Element-wise scaling by a complex scalar.
+   * @param scalar Factor applied to every matrix entry.
+   * @return Scaled copy of this matrix.
+   */
+  [[nodiscard]] DynamicMatrix operator*(const Complex& scalar) const;
+
+  /**
+   * @brief Element-wise in-place scaling by a complex scalar.
+   * @param scalar Factor applied to every matrix entry.
+   * @return Reference to this matrix.
+   */
+  DynamicMatrix& operator*=(const Complex& scalar);
+
+  /**
+   * @brief Checks whether this matrix is approximately the identity.
+   * @param tol Maximum allowed complex modulus of each off-diagonal entry and
+   * each diagonal deviation from one.
+   * @return True when the matrix is close to the identity.
+   */
+  [[nodiscard]] bool isIdentity(double tol = MATRIX_TOLERANCE) const;
+
 private:
   struct Impl;
   std::unique_ptr<Impl> impl_;
@@ -558,4 +768,31 @@ inline constexpr bool
     std::disjunction_v<std::is_same<T, Matrix1x1>, std::is_same<T, Matrix2x2>,
                        std::is_same<T, Matrix4x4>,
                        std::is_same<T, DynamicMatrix>>;
+
+/// Scalar-on-the-left multiply `scalar * matrix` (commutes with the member
+/// `matrix * scalar`). Provided so generic code can scale a matrix from
+/// either side.
+[[nodiscard]] Matrix2x2 operator*(const Complex& scalar,
+                                  const Matrix2x2& matrix);
+/// @copydoc operator*(const Complex&, const Matrix2x2&)
+[[nodiscard]] Matrix4x4 operator*(const Complex& scalar,
+                                  const Matrix4x4& matrix);
+/// @copydoc operator*(const Complex&, const Matrix2x2&)
+[[nodiscard]] DynamicMatrix operator*(const Complex& scalar,
+                                      const DynamicMatrix& matrix);
+
+/**
+ * @brief Eigenvalues and eigenvectors of a real symmetric `4x4` matrix.
+ *
+ * `eigenvalues` are sorted ascending and `eigenvectors` holds the
+ * corresponding orthonormal eigenvectors as columns (column `j` is the
+ * eigenvector for `eigenvalues[j]`).
+ */
+struct SymmetricEigen4 {
+  /// Eigenvalues in ascending order.
+  std::array<double, 4> eigenvalues{};
+  /// Orthonormal eigenvectors as columns (column `j` matches `eigenvalues[j]`).
+  Matrix4x4 eigenvectors{};
+};
+
 } // namespace mlir::qco