Phys 142 - A Complex Matrix Multiplication Program

Note: My code here starts with:

#include <iostream>
#include <vector>
#include <complex>
using namespace std;

A good introduction for iostream is cplusplus.com’s Basic Input/Output article. It is basically printf, except you write cout<<"Number is: "<<number<<endl; to output something instead of printf("Decimals: %d", 1000);.

A quick introduction to vector is on cprogramming.com. It means you can create dynamically resizable arrays and return arrays without worrying about those pesky new and delete keywords.

A quick reference for complex is the cppreference reference page (code examples at bottom).

Finally: using namespace std means we can write cout<<"Hello"<<endl; instead of std::cout<<"Hello"<<std::endl;

If you object to my use of std::complex or row-major matrices, I am doing it this way so that things play nice with BLAS and CUDA!

Intro

I want to come up with an instructive example of raising some approximately unitary matrix to a large power.

Let’s work with the matrix:

We’ll take $\theta=2\pi$. This matrix, raised to the $N$th power (with $\theta=2\pi$), is approximately the negative identity matrix. (This is exact as $N\to \infty$). It’s approximately unitary because it’s of the form $M=1+i\varepsilon A$ with $A$ Hermitian, so that $M M^\dagger=(1+i\varepsilon A)(1-i \varepsilon A)=1+O(\varepsilon^2)$. Let me first give the whole commented code snippet, then explain it. You can compile it by typing g++ -o matrix matrix.cpp and run it by typing ./matrix in your terminal.

#include <iostream>
#include <complex>
#include <vector>
using namespace std;

//working with 2 by 2 matrices
const int D=2;

//This line lets us write "cdouble x;" instead of "complex<double> x",
//and "vector<cdouble>" arr; instead of "vector<complex<double> > arr;"
typedef complex<double> cdouble;

//With this function we can write "A=matrixMultiply(B,C);". 
vector<cdouble> matrixMultiply(vector<cdouble> m1,vector<cdouble> m2);

int main() {
    //Declare a complex datatype
    cdouble I(0,1);
    double theta=2*3.1415926;
    int N=1000;

    //Construct the matrix we'll be multiplying by over and over again. 
    vector<cdouble> m0;
    m0.push_back(1);    
    m0.push_back(I*(theta/(2*N)));  
    m0.push_back(I*(theta/(2*N)));  
    m0.push_back(1);    
    
    //Construct the identity matrix.
    vector<cdouble> matrix;
    matrix.push_back(1);
    matrix.push_back(0);
    matrix.push_back(0);
    matrix.push_back(1);

    //left-multiply by the matrix N times.
    for(int n=0;n<N;n++){
        matrix=matrixMultiply(m0,matrix);
    }
    //Print the matrix in a readable format.
    cout<<"["<<matrix[0]<<"\t"<<matrix[1]<<"]"<<endl;
    cout<<"["<<matrix[2]<<"\t"<<matrix[3]<<"]"<<endl;
    return 0;
}

vector<cdouble> matrixMultiply(vector<cdouble> m1,vector<cdouble> m2){
    vector<cdouble> matrixnew(D*D,0);
    //loop over the rows and columns of matrixnew
    for(int i=0;i<D;i++){
        for(int j=0;j<D;j++){
            //zero the i,jth element and then sum over k
            matrixnew[D*i+j]=0;
            for(int k=0;k<D;k++){
                //repeated index summation convention:
                //mnew_ij=sum_k(m_ik*m_kj)
                matrixnew[D*i+j]+=m1[D*i+k]*m2[D*k+j];
            }
        }
    }   
    //return the vector object.
    return matrixnew;
}

Matrix multiplication and row-major matrices

When we implement this using BLAS, the easiest thing to do is put the matrix in row-major format. This means instead of a two-dimensional array (of type double** or of type vector< vector < double> >), we have a one-dimensional array (of type double*, or vector<double>). The way to do this is to concatenate all the rows of $M_0$ together. The first row of $M_0$ is $\begin{bmatrix} 1 & -\frac{i \theta}{2N} \end{bmatrix} $. The second row of $M_0$ is $\begin{bmatrix} -\frac{i \theta}{2N} & 1 \end{bmatrix} $. The matrix is represented in row-major format as $\begin{bmatrix} 1 & -\frac{i \theta}{2N} &-\frac{i \theta}{2N} & 1\end{bmatrix} $.

This is, for example, how the bitmap datatype works. You have a raster image $W$ pixels wide and $H$ pixels tall, stored in a string of bits (AKA a one-dimensional array). There are $W\cdot H$ pixels, and n ranges from 0 to $WH-1$ inclusive. The first row of pixels are pixels 0<=n<W. The second row has W<=n<2*W, and so on. In general, the pixel at row i, column j, is at location W*i+j. So this means we can access the ith row and jth column of a matrix mat by typing in mat[W*i+j].

We want to compute the matrix product of square $D$-dimensional matrices, $C=A * B$. Using index notation, we have $C_{ij}=\sum_k A_{ik} B_{kj}$ If A and B are std::vector objects in row-major format, this becomes a sum over k of A[D*i+k]*B[D*k+j].

In an algorithm: we do for loops over i and j, and then sum over the k variable. Easy peasy:

 
vector<double> matrixMultiply(vector<double> m1,vector<double> m2){
    vector<double> matrixnew(D*D,0);
    for(int i=0;i<D;i++){
        for(int j=0;j<D;j++){
            matrixnew[D*i+j]=0;
            for(int k=0;k<D;k++){
                //repeated index summation convention:
                //mnew_ij=sum_k(m_ik*m_kj)
                matrixnew[D*i+j]+=m1[D*i+k]*m2[D*k+j];
            }
        }
    }   
    return matrixnew;
}

Complex numbers and complex matrices

The complex number type behaves about how you’d expect it to. You can have floating point or integer coefficient complex numbers, but it’s easiest to focus on double floating point numbers. Here is an example program showing the capabilities of complex numbers. You can compile it on the lab computers using g++ -o complex complex.cpp, and run it with ./complex.

 
#include <iostream>
#include <complex>
using namespace std;

int main() {
    complex<double> x(0.7071,0.7071);
    complex<double> I(0,1);

    cout<<"0.707+0.707*i squared is: "<<(x*x)<<endl;
    cout<<"i squared is: "<<(I*I)<<endl;
    cout<<"(i+1)/(0.707+0.707*i) is: "<<((I+1.0)/x)<<endl;
    return 0;
}

Note: complex<double> can be multiplied and divided by type double. It doesn’t work so well with type int! Hence why I wrote I+1.0 instead of I+1.

Typing in complex<double> gets old fast, so we can create a typedef shortcut:

 
typedef complex<double> cdouble;

int main() {
    cdouble x(0.7071,0.7071);
    cout<<"0.707+0.707*i squared is: "<<(x*x)<<endl;
    return 0;
}

Finally, let’s construct that array from above using std::vector, storing it in row-major format:

Our construction is then:

#include <iostream>
#include <complex>
#include <vector>
using namespace std;

typedef complex<double> cdouble;

int main() {
    cdouble I(0,1);
    double theta=2*3.1415926;
    int N=100;
    
    vector<cdouble> m0;
    m0.push_back(1); //element 0
    m0.push_back(-I*(theta/(2*N))); //element 1
    m0.push_back(-I*(theta/(2*N))); //element 2
    m0.push_back(1); //element 3

    //print out the element in the first row and second column.
    cout<<m0[1]<<endl;

    return 0;
}

Where the matrix comes from

The example we used was a spin one-half particle. From introductory quantum mechanics, the matrices $S^k=\frac{i}{2}\sigma^k$ are a spin one-half representation of the su(2) algebra $[S^i,S^j]=i\varepsilon _ {ink}S^k$. Spin one-half is weird: objects turn into minus themselves after a $360^\circ$ rotation. Let’s demonstrate this fact numerically.

Say we have a spin one-half particle in state $|\psi\rangle$. A rotation of this particle about the x axis by angle theta can be enacted by calculating $e^{-i\theta S^x}|\psi\rangle$. Expanding…

Sigma x is the matrix:

So the matrix we want to raise to the nth power is