Suppose that a matrix 



M


{\displaystyle M}

 contains complex numbers. Then the conjugate transpose of 



M


{\displaystyle M}

 is a matrix 




M

H




{\displaystyle M^{H}}

 containing the complex conjugates of the matrix transposition of 



M


{\displaystyle M}

.

 



(

M

H



)

j
i


=



M

i
j


¯




{\displaystyle (M^{H})_{ji}={\overline {M_{ij}}}}




This means that row 



j


{\displaystyle j}

, column 



i


{\displaystyle i}

 of the conjugate transpose equals the

complex conjugate of row 



i


{\displaystyle i}

, column 



j


{\displaystyle j}

 of the original matrix.



In the next list, 



M


{\displaystyle M}

 must also be a square matrix.

 A Hermitian matrix equals its own conjugate transpose: 




M

H


=
M


{\displaystyle M^{H}=M}

.
 A normal matrix is commutative in multiplication with its conjugate transpose: 




M

H


M
=
M

M

H




{\displaystyle M^{H}M=MM^{H}}

.
 A unitary matrix has its inverse equal to its conjugate transpose: 




M

H


=

M

−
1




{\displaystyle M^{H}=M^{-1}}

. 
 This is true iff 




M

H


M
=

I

n




{\displaystyle M^{H}M=I_{n}}

 and iff 



M

M

H


=

I

n




{\displaystyle MM^{H}=I_{n}}

, where 




I

n




{\displaystyle I_{n}}

 is the identity matrix.


Task

Given some matrix of complex numbers, find its conjugate transpose.

Also determine if the matrix is a:

 Hermitian matrix,
 normal matrix, or 
 unitary matrix.


See also

 MathWorld entry: conjugate transpose
 MathWorld entry: Hermitian matrix
 MathWorld entry: normal matrix
 MathWorld entry: unitary matrix


#include <cassert>
#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <vector>
 
template <typename scalar_type> class complex_matrix {
public:
    using element_type = std::complex<scalar_type>;
 
    complex_matrix(size_t rows, size_t columns)
        : rows_(rows), columns_(columns), elements_(rows * columns) {}
 
    complex_matrix(size_t rows, size_t columns, element_type value)
        : rows_(rows), columns_(columns), elements_(rows * columns, value) {}
 
    complex_matrix(size_t rows, size_t columns,
        const std::initializer_list<std::initializer_list<element_type>>& values)
        : rows_(rows), columns_(columns), elements_(rows * columns) {
        assert(values.size() <= rows_);
        size_t i = 0;
        for (const auto& row : values) {
            assert(row.size() <= columns_);
            std::copy(begin(row), end(row), &elements_[i]);
            i += columns_;
        }
    }
 
    size_t rows() const { return rows_; }
    size_t columns() const { return columns_; }
 
    const element_type& operator()(size_t row, size_t column) const {
        assert(row < rows_);
        assert(column < columns_);
        return elements_[row * columns_ + column];
    }
    element_type& operator()(size_t row, size_t column) {
        assert(row < rows_);
        assert(column < columns_);
        return elements_[row * columns_ + column];
    }
 
    friend bool operator==(const complex_matrix& a, const complex_matrix& b) {
        return a.rows_ == b.rows_ && a.columns_ == b.columns_ &&
               a.elements_ == b.elements_;
    }
 
private:
    size_t rows_;
    size_t columns_;
    std::vector<element_type> elements_;
};
 
template <typename scalar_type>
complex_matrix<scalar_type> product(const complex_matrix<scalar_type>& a,
                                    const complex_matrix<scalar_type>& b) {
    assert(a.columns() == b.rows());
    size_t arows = a.rows();
    size_t bcolumns = b.columns();
    size_t n = a.columns();
    complex_matrix<scalar_type> c(arows, bcolumns);
    for (size_t i = 0; i < arows; ++i) {
        for (size_t j = 0; j < n; ++j) {
            for (size_t k = 0; k < bcolumns; ++k)
                c(i, k) += a(i, j) * b(j, k);
        }
    }
    return c;
}
 
template <typename scalar_type>
complex_matrix<scalar_type>
conjugate_transpose(const complex_matrix<scalar_type>& a) {
    size_t rows = a.rows(), columns = a.columns();
    complex_matrix<scalar_type> b(columns, rows);
    for (size_t i = 0; i < columns; i++) {
        for (size_t j = 0; j < rows; j++) {
            b(i, j) = std::conj(a(j, i));
        }
    }
    return b;
}
 
template <typename scalar_type>
std::string to_string(const std::complex<scalar_type>& c) {
    std::ostringstream out;
    const int precision = 6;
    out << std::fixed << std::setprecision(precision);
    out << std::setw(precision + 3) << c.real();
    if (c.imag() > 0)
        out << " + " << std::setw(precision + 2) << c.imag() << 'i';
    else if (c.imag() == 0)
        out << " + " << std::setw(precision + 2) << 0.0 << 'i';
    else
        out << " - " << std::setw(precision + 2) << -c.imag() << 'i';
    return out.str();
}
 
template <typename scalar_type>
void print(std::ostream& out, const complex_matrix<scalar_type>& a) {
    size_t rows = a.rows(), columns = a.columns();
    for (size_t row = 0; row < rows; ++row) {
        for (size_t column = 0; column < columns; ++column) {
            if (column > 0)
                out << ' ';
            out << to_string(a(row, column));
        }
        out << '\n';
    }
}
 
template <typename scalar_type>
bool is_hermitian_matrix(const complex_matrix<scalar_type>& matrix) {
    if (matrix.rows() != matrix.columns())
        return false;
    return matrix == conjugate_transpose(matrix);
}
 
template <typename scalar_type>
bool is_normal_matrix(const complex_matrix<scalar_type>& matrix) {
    if (matrix.rows() != matrix.columns())
        return false;
    auto c = conjugate_transpose(matrix);
    return product(c, matrix) == product(matrix, c);
}
 
bool is_equal(const std::complex<double>& a, double b) {
    constexpr double e = 1e-15;
    return std::abs(a.imag()) < e && std::abs(a.real() - b) < e;
}
 
template <typename scalar_type>
bool is_identity_matrix(const complex_matrix<scalar_type>& matrix) {
    if (matrix.rows() != matrix.columns())
        return false;
    size_t rows = matrix.rows();
    for (size_t i = 0; i < rows; ++i) {
        for (size_t j = 0; j < rows; ++j) {
            if (!is_equal(matrix(i, j), scalar_type(i == j ? 1 : 0)))
                return false;
        }
    }
    return true;
}
 
template <typename scalar_type>
bool is_unitary_matrix(const complex_matrix<scalar_type>& matrix) {
    if (matrix.rows() != matrix.columns())
        return false;
    auto c = conjugate_transpose(matrix);
    auto p = product(c, matrix);
    return is_identity_matrix(p) && p == product(matrix, c);
}
 
template <typename scalar_type>
void test(const complex_matrix<scalar_type>& matrix) {
    std::cout << "Matrix:\n";
    print(std::cout, matrix);
    std::cout << "Conjugate transpose:\n";
    print(std::cout, conjugate_transpose(matrix));
    std::cout << std::boolalpha;
    std::cout << "Hermitian: " << is_hermitian_matrix(matrix) << '\n';
    std::cout << "Normal: " << is_normal_matrix(matrix) << '\n';
    std::cout << "Unitary: " << is_unitary_matrix(matrix) << '\n';
}
 
int main() {
    using matrix = complex_matrix<double>;
 
    matrix matrix1(3, 3, {{{2, 0}, {2, 1}, {4, 0}},
                          {{2, -1}, {3, 0}, {0, 1}},
                          {{4, 0}, {0, -1}, {1, 0}}});
 
    double n = std::sqrt(0.5);
    matrix matrix2(3, 3, {{{n, 0}, {n, 0}, {0, 0}},
                          {{0, -n}, {0, n}, {0, 0}},
                          {{0, 0}, {0, 0}, {0, 1}}});
 
    matrix matrix3(3, 3, {{{2, 2}, {3, 1}, {-3, 5}},
                          {{2, -1}, {4, 1}, {0, 0}},
                          {{7, -5}, {1, -4}, {1, 0}}});
 
    test(matrix1);
    std::cout << '\n';
    test(matrix2);
    std::cout << '\n';
    test(matrix3);
    return 0;
}