Linear algebra matrix calculations

determinant.rb

require "matrix"
require "benchmark"

def mprint(m)
#pretty-print matrices
  puts m.map{|z| z.map{|q| q.round(4).to_s.rjust(7," ")}.join(", ")}.join("\n")
end


def minor(m,i,j)
#returns minor matrix (comatrix)
  r=m.size
  c=m[0].size
  ret=Array.new(r-1){Array.new(c-1){nil}}
  
  for n in 0..i-1 do
    for k in 0..j-1 do
      ret[n][k]=m[n][k]
    end
  end
  
  for n in 0..i-1 do
    for k in j+1..c-1 do
      ret[n][k-1]=m[n][k]
    end
  end
  
    
  for n in i+1..r-1 do
    for k in 0..j-1 do
      ret[n-1][k]=m[n][k]
    end
  end
  
  for n in i+1..r-1 do
    for k in j+1..c-1 do
      ret[n-1][k-1]=m[n][k]
    end
  end
 return ret
end

def recdet(m)
#recursively compute determinant by Laplace cofactor expansion, slow
  r=m.size
  c=m[0].size
  if r==c && r==2 then
    return m[0][0]*m[1][1]-m[0][1]*m[1][0]
  else
  
    i=rand(r)
    d=m[i].collect.with_index.map{|z,j| (((-1)**(i+j))*z)*det(minor(m,i,j))}.inject(:+)
	    
  end
  
  return d  
end

def det(m)
#compute determinant by converting matrix to low triangle and then taking the the diagonal product, faster
  ret=triangle(m)
  return ((-1)**ret[1])*diagprod(ret[0])
end

def swap(m,i,j)
#swaps i-th and j-th rows in the matrix
  ret=m.map(&:clone)
  temp=ret[i]
  ret[i]=ret[j]
  ret[j]=temp
  return ret
end

def mul(m,r,a)
#multiply one row with a value
  ret=m.map(&:clone)
  ret[r]=m[r].map{|z| z*a}
  return ret
end


def addmul(m,r1,r2,a)
#adds multiplied row r2 by value a to r1
   ret=m.map(&:clone)
   c=m[0].size
   for i in 0..c-1 do
     ret[r1][i]=ret[r1][i]+ret[r2][i]*a
   end
   return ret
end

def triangle(m)
#makes matrix low triangle form, row spaws recorded
  ret=m.map(&:clone)
  r=ret.size
  swaps=0
  for j in 0..r-2
    for i in j+1..r-1 do
      maxind=ret.collect.with_index.to_a[j..-1].max_by{|z,w| z[j].abs}[1]
      if j!=maxind then
        ret=swap(ret,j,maxind)
        swaps+=1
      end
      ret=addmul(ret,i,j,-ret[i][j]/ret[j][j])
    end
  end
  
  return [ret,swaps]
end

def inverse(m)
#calculate inverse matrix by Gauss elimination and partial pivoting
  ret=m.map(&:clone)
  r=ret.size
  
  ones=Array.new(r){Array.new(r){0.0}}
  
  for i in 0..r-1
    ones[i][i]=1.0
  end
    
  for j in 0..r-2
    for i in j+1..r-1 do
      maxind=ret.collect.with_index.to_a[j..-1].max_by{|z,w| z[j].abs}[1]
      if j!=maxind then
        ret=swap(ret,j,maxind)
        ones=swap(ones,j,maxind)
        
      end
      f=-ret[i][j].fdiv(ret[j][j])
      ret=addmul(ret,i,j,f)
      ones=addmul(ones,i,j,f)
      
    end
  end
  
  for i in 0..r-2
    f=1.0/ret[i][i]
    ret=mul(ret,i,f)
    ones=mul(ones,i,f)
    
    for j in i+1..r-1 do
      
      maxind=ret.collect.with_index.to_a[i..-1].max_by{|z,w| z[i].abs}[1]
      if i!=maxind then
        ret=swap(ret,i,maxind)
        ones=swap(ones,i,maxind)
        
      end
      f=-ret[i][j].fdiv(ret[j][j])
      ret=addmul(ret,i,j,f)
      ones=addmul(ones,i,j,f)
      
    end
    
  end
  
  
  f=1.0/ret[r-1][r-1]
  ret=mul(ret,r-1,f)
  ones=mul(ones,r-1,f)
  
  
  return ones
end



def diagprod(m)
#calculate the diagonal product of a low triangle matrix
  r=m.size
  prod=1.0
  for i in 0..r-1
    prod*=m[i][i]
  end
  return prod
end

def mmul(m1,m2)
#multiply two matrices
  r=m1.size
  c=m2[0].size
  if r!=c then
    return nil
  end
  temp=m2.transpose
  ret=Array.new(r){Array.new(c){nil}}
  for i in 0..r-1 do
    for j in 0..c-1 do
      ret[i][j]=m1[i].zip(temp[j]).map{|z| z[0]*z[1]}.inject(:+)
    end
  end
  
  return ret
end

def randm(n)
#generate random matrix
  m=Array.new(n){Array.new(n){(rand).round(4)}}
end



n=8

t=Benchmark.measure do 
  for i in 0..1000 do
    m=Array.new(n){Array.new(n){1-2*rand}}
    m2=Matrix[*m]
    print m2.det 
    print " ~ "
  end
end

t2=Benchmark.measure do 
  for i in 0..1000 do
    m=Array.new(n){Array.new(n){1-2*rand}}
    m2=Matrix[*m]
    mprint m
    print det(m)
    print " ~ "
  end
end

puts 
puts t
puts t2