Search code examples
rmatrixeigenvector

Find Matrix which has orthogonal eigenvectors


I would like to find an n x n matrix S, whose eigenvectors are orthogonal which also has non-repeated integer eigenvalues. A test for an orthogonal matrix Q is QQ^T=I. I have written a function to test 1e5 random matrices. How can I return the matrix if it fulfils the criteria, but throw it away if it does not?

1. For a single test matrix:

require(matrix)

s<-matrix(sample(1:20,9,replace=TRUE),ncol=3) #Random Test matrix
s[lower.tri(s)] = t(s)[lower.tri(s)] #Make it symmetrical

eigen(s)->dC #Find eigenvectors and eigenvalues of S.

#Test orthogonality of matrix of eigenvectors (Q) by QQ^T=Identity

sum(dC$vectors%*%t(dC$vectors)-diag(1,3))<1e-6 #Gives TRUE or FALSE

#Test for integer eigenvalues:

   if(dC$values%%1==0){return(TRUE)} #Want integer eigenvalues
      else {return(FALSE)}

2. Trying to test many matrices

#Create some symmetrical matrices 
testMatrices<-list() #List of matrices
for (k in 1:1e5) {
  (s<-matrix(sample(1:20,9,replace=TRUE),ncol=3))
 s[lower.tri(s)] = t(s)[lower.tri(s)] #Create triangular matrix
  testMatrices[[k]] <- s   
};rm(k)

eigenDecompFunction<-function(s){
  eigen(s)->dC

  if(dC$values%%1==0){return(TRUE)} #Want integer eigenvalues
  else {return(FALSE)}

  
 #Test for orthogonal matrix QQ^T=I
  if (sum(dC$vectors%*%t(dC$vectors)-diag(1,3))<1e-6){
    return(s)
  }else {return(FALSE)}
}

Solution

  • It is much more efficient to generate an 'easy' matrix with the desired eigenvectors and then apply random rotations. Using these two Mathematics Stack Exchange answers to brush up my rusty linear algebra:

    n <- 5
    v <- 1:n        ## specify eigenvalues (one simple example!)
    m <- diag(1:n)  ## create a diag matrix with these elements
    set.seed(101)
    ## from Math SE link #1: could also use `pracma::randortho(n, type = "unitary")
    r1 <- matrix(rnorm(n^2), n, n)
    q1 <- qr(r1)
    ## from Math SE link #2: this is t(Q) %*% m %*% Q
    m2 <- qr.qty(q1, qr.qy(q1, m))
    

    Now checking with eigen(m2)$values, the result is 5 4 3 2 1 (original eigenvalues, sorted in decreasing order as is standard for eigenvalues)


    The solution above gives random real-valued (or rather floating-point-valued) matrices with integral eigenvalues (or rather floating-point values that might be integer-valued or very close). If you want an integer-valued matrix with integer eigenvalues, I have some bad news for you:

    Our main result affirms and strengthens a conjecture made in [4]: for any n > 2, the probability that a random n x n integer matrix has even a single integer eigenvalue is 0.

    Martin, Greg, and Erick B. Wong. “Almost All Integer Matrices Have No Integer Eigenvalues.” The American Mathematical Monthly 116, no. 7 (2009): 588–97.

    And then I have some worse news, which is that the rational eigenvalues of an integer-valued matrix are integral. In combination with the previous result, that means that whether you want integer-valued or just rational eigenvalues, you have a mathematically impossible problem to solve. That doesn't mean there aren't special forms of integer matrices that will satisfy your constraints (potentially providing infinitely many solutions) — but unless I'm missing something it means that your strategy of sampling random integer matrices and hoping that some of them will have rational or integer eigenvalues is effectively doomed ...