Class: AppMath::Mat

Inherits:

Object

• Object
• AppMath::Mat

Includes:

Comparable, Enumerable

Defined in:

lib/linalg.rb

Overview

Matrix space of arbitrary dimension. The intended usage is that the elements of a matrix are all either real or complex. Since one is allowed to change any matrix element into any object there is no guaranty for type-uniformity of the elements of a matrix.

Instance Attribute Summary

• #x ⇒ Object

  Returns the value of attribute x.

Class Method Summary

• .svdcmp(a, w, v) ⇒ Object

  Singular value decomposition.

• .test(n0, verbose = true, complex = false) ⇒ Object

  Consistency test of class Mat.

• .tob(n, i, complex = false) ⇒ Object

  Generates a test object, here a n times n matrix with random elements.

Instance Method Summary

• #*(v) ⇒ Object

  Multiplication of a Mat with either a Mat, Vec, or Numeric.

• #+(v) ⇒ Object

  Returns self + v, where v is a Mat.

• #-(v) ⇒ Object

  Returns self - v , where v is a Mat.

• #-@ ⇒ Object

  Unary minus operator.

• #<=>(v) ⇒ Object

  The order relation is here lexicographic ordering of lists.

• #[](i) ⇒ Object

  Reading row vectors.

• #[]=(i, a) ⇒ Object

  Setting row vectors.

• #abs ⇒ Object

  Absolute value, always real.

• #abs2 ⇒ Object

  Square of absolute value.

• #clone ⇒ Object

  Returns an independent copy of self.

• #complex? ⇒ Boolean
• #conj ⇒ Object

  Returns the Hermitian conjugate of self.

• #dim ⇒ Object

  Returns the ‘dimension’ of the matrix, i.e.

• #dim1 ⇒ Object

  Let self be a (m,n)-matrix (also called a m times n matrix) then dim1 == m.

• #dim2 ⇒ Object

  Let self be a (m,n)-matrix (also called a m times n matrix) then dim2 == n.

• #dis(v) ⇒ Object

  Relative distance between matrices.

• #each ⇒ Object
• #initialize(*arg) ⇒ Mat constructor

  These are the 5 mehods to generate a matrix via ‘new’ m1 = Mat.new(aMat) m2 = Mat.new(anArrayOfVec) m3 = Mat.new(aVec) m4 = Mat.new(aPositiveInteger, aRealOrComplex) m5 = Mat.new(aPositiveInteger, aPositiveInteger, aRealOrComplex) Here, m1 is a copy of aMat, m2 is a matrix which has as row vectors, the components of anArrayOfVec.

• #inv ⇒ Object

  In most cases (‘up to a subst of measure zero’) the pseudo-inverse is also the inverse.

• #prn(name) ⇒ Object

  Prints the content of self and naming the output.

• #pseudo_inv(acc = 0) ⇒ Object

  Returns the pseudo-inverse (als known as Penrose inverse) of self.

• #s!(i, j, a) ⇒ Object

  The ‘s’ stands for ‘set (value)’ and the ‘!’ this is a method by which self changes (non-constant or mutating method).

• #spr(v) ⇒ Object

  Scalar product of matrices.

• #square? ⇒ Boolean
• #to_0 ⇒ Object

  ‘to zero’ Returns a matrix with he same dimensions as self, but with all matrix elements set to zero.

• #to_1 ⇒ Object

  ‘to one’.

• #to_s ⇒ Object

  Returns a string which consists of a list of the strings which represent the row-vectors.

• #trp ⇒ Object

  Returns the transposed of self.

Constructor Details

#initialize(*arg) ⇒ Mat

These are the 5 mehods to generate a matrix via ‘new’

m1 = Mat.new(aMat)
m2 = Mat.new(anArrayOfVec)
m3 = Mat.new(aVec)
m4 = Mat.new(aPositiveInteger, aRealOrComplex)
m5 = Mat.new(aPositiveInteger, aPositiveInteger, aRealOrComplex)

Here, m1 is a copy of aMat, m2 is a matrix which has as row vectors, the components of anArrayOfVec. If these vectors have not all he same dimension, failure results; m3 is a square matrix in which only the main diagonal may have non-zero elements, and in which ths diagonal is given as aVec; m4 is a square matrix with the dimension given by the first argument, and with all matrix elements equal to the second argment; m5 is a rectangular matrix with dim1 and dim2 given by the first and the second argument, and with all matrix elements equal to the third argument.

# File 'lib/linalg.rb', line 429

def initialize(*arg)
  case arg.size
  when 0
    @x = Array.new
  when 1 # ok if this is a Matrix or an array of Vectors
    a0 = arg[0]
    if a0.is_a?(Mat)
      @x = Array.new(a0.x)
    elsif a0.is_a?(Array)
      n = a0.size
      if n.zero?
        @x = Array.new
      else
        misfit = 0
        a0.each{|c|
          misfit += 1 unless c.is_a?(Vec)
        }
        fail "input must consist of Vec-objects" unless misfit.zero?
        misfit2 = 0
        d2 = a0[1].dim
        a0.each{|c|
          misfit2 += 1 unless c.dim == d2
        }
        fail "input Vec-objects must agree in dimension" unless 
          misfit.zero?
        @x = a0.clone
      end
    elsif a0.is_a?(Vec) # make a diagonal matrix
      n = a0.dim
      if n.zero?
        @x = Array.new
      else
        c = a0[1].to_0
        vc = Vec.new(n,c)
        @x = Array.new(n,vc)
        for i in 1..n
          s!(i,i,a0[i])
        end
      end
    else
      fail "no reasonable construction available for this argument"
    end
  when 2 # make a square matrix, the diagonal filled with one element
    # (all others zero)
    n = arg[0]
    a = arg[1]
    zero = a.to_0
    vc = Vec.new(n,zero)
    @x = Array.new(n,vc)
    for i in 1..n
      vi = Vec.new(vc)
      vi[i] = a
      @x[i-1] = vi
    end
  when 3 # make rectangular matrix filled with one element
    n1 = arg[0]
    fail "first argument must be integer" unless n1.integer?
    n2 = arg[1]
    fail "second argument must be integer" unless n2.integer?
    a = arg[2]
    vc = Vec.new(n2,a)
    @x = Array.new(n1,vc)
  else
    fail "no construction for more than 3 arguments"
  end
end

Instance Attribute Details

#x ⇒ Object

Returns the value of attribute x.

# File 'lib/linalg.rb', line 396

def x
  @x
end

Class Method Details

.svdcmp(a, w, v) ⇒ Object

Singular value decomposition. Slightly modified fom Press et al. Only needed as a algorithmic tool. The method for the end-user is method pseudo_inv.

# File 'lib/linalg.rb', line 538

def Mat.svdcmp(a, w, v)

  m = a.dim1; n = a.dim2
  fail "svdcmp: bad frame of a" unless m >= n
  fail "svdcmp: bad frame of w" unless n == w.dim
  fail "svdcmp: bad frame of v" unless v.dim1 == n && v.dim2 == n
  fail "svdcmp: dim = 0 as input" if m.zero? || n.zero?

  iter_max=40
  a11 = a[1][1]
  zero =a11.to_0
  one = a11.to_1
  two = one + one
  rv1 = Vec.new(n,zero)
  g = zero; scale = zero; anorm = zero
  for i in 1..n 
    l = i + 1
    rv1[i] = scale * g
    g = zero; s = zero; scale = zero
    if i <= m
      for k in i..m; scale += a[k][i].abs; end
      if scale.nonzero?
        for k in i..m
          aki = a[k][i]
          aki /= scale
          a.s!(k,i,aki)
          s += aki * aki
        end
        f = a[i][i]
        g = - Basics.sign2(s.sqrt,f)
        h = f * g - s
        a.s!(i,i,f - g)
        for j in l..n
          s = zero
          for k in i..m; s += a[k][i] * a[k][j]; end
          f = s/h
          for  k in i..m 
            akj = a[k][j]
            akj += f * a[k][i]
            a.s!(k,j,akj)
          end
        end # for j in l..n
        for k in i..m
          aki = a[k][i]
          aki *= scale
          a.s!(k,i,aki)
        end
      end # scale != zero
    end # i <= m
    w[i] = scale * g
    g = zero; s = R.c0; scale = zero
    if i <= m && i != n
      for k in l..n; scale += a[i][k].abs; end
      if scale.nonzero?
        for k in l..n
          aik = a[i][k]
          aik /= scale
          a.s!(i,k,aik)
          s += aik * aik
        end
        f = a[i][l]
        g = - Basics.sign2(s.sqrt,f)
        h = f * g - s
        a.s!(i,l,f - g)
        for k in l..n; rv1[k] = a[i][k]/h ; end
        for j in l..m
          s = zero
          for k in l..n; s += a[j][k] * a[i][k]; end
          for k in l..n
            ajk = a[j][k]
            ajk += s * rv1[k]
            a.s!(j,k,ajk)
          end
        end # for j in l..m
        for k in l..n; aik = a[i][k]; aik *= scale; a.s!(i,k,aik); end
      end # if scale != zero
    end # if i <= m && i != n
    anorm = Basics.sup(anorm,w[i].abs + rv1[i].abs)
  end # for i in 1..n 
  i = n
  while i >= 1
    if i < n
      if g.nonzero?
        for j in l..n; v.s!(j,i, (a[i][j]/a[i][l])/g); end
        for j in l..n
          s = zero
          for k in l..n; s += a[i][k] * v[k][j]; end
          for k in l..n
            vkj =v[k][j]
            vkj += s * v[k][i]
            v.s!(k,j,vkj)
          end
        end # for j in l..n
      end # if g.notzero!
      for j in l..n; v.s!(i,j,zero); v.s!(j,i,zero); end
    end # if i< n
    v.s!(i,i,one)
    g = rv1[i]
    l = i
    i -= 1
  end # while i >= 1

  i = Basics.inf(m,n)
  while i >= 1
    l = i + 1
    g = w[i]
    for j in l..n; a.s!(i,j,zero); end
    if g.nonzero?
      g = one/g
      for j in l..n
        s = zero
        for k in l..m; s += a[k][i] * a[k][j]; end
        f = (s/a[i][i]) * g
        for k in i..m
          akj = a[k][j]; akj += f * a[k][i]; a.s!(k,j,akj)
        end
      end # for j in l..n
      for j in i..m; aji = a[j][i]; aji *= g; a.s!(j,i,aji); end
    else
      for j in i..m; a.s!(j,i,zero); end
    end # if g.nonzero?  
    aii = a[i][i]; aii += one; a.s!(i,i,aii)
    i -= 1
  end # while i >= 1

  k = n
  while k >=1
    for its in 1..iter_max
      flag = 1
      l = k
      while l >= 1
        nm = l - 1
        if rv1[l].abs + anorm == anorm
          flag = 0
          break
        end # if rv1[l].abs + anorm == anorm
        break if w[nm].abs + anorm == anorm
        l -= 1
      end #  while l >= 1
      if flag.nonzero?
        c = zero
        s = one
        for i in l..k
          f = s * rv1[i]
          rv1[i] = c * rv1[i]
          break if f.abs + anorm == anorm
          g = w[i]
          h = f.hypot(g)
          w[i] = h
          h = one/h
          c = g * h
          s = -f * h
          for j in 1..m
            y = a[j][nm]; z = a[j][i]; 
            a.s!(j,nm,y*c+z*s); a.s!(j,i,z*c-y*s)
          end # for j in 1..m
        end # for i in l..k
      end # if flag.nonzero?
      z = w[k]
      if l == k
        if z < zero
          w[k] = -z
          for j in 1..n; v.s!(j,k,-v[j][k]); end 
        end # if z < zero
        break
      end # if l == k
      fail "no convergence in #{iter_max} iterations" if its == iter_max
      x = w[l]; nm = k - 1; y = w[nm]; g = rv1[nm]; h = rv1[k]
      f = ((y-z) * (y+z) + (g-h) * (g+h))/(h*y*two)
      g = f.hypot(one)
      f = ((x-z)*(x+z)+h*((y/(f+Basics.sign2(g,f)))-h))/x;
      c = one; s = one
      for j in l..nm
        i = j + 1; g = rv1[i]; y =w[i]; h = s * g; g = c * g
        z = f.hypot(h)
        rv1[j] = z
        c = f/z
        s = h/z
        f = x*c+g*s; g = g*c-x*s; h=y*s; y *= c;
        for jj in 1..n
          x=v[jj][j]; z=v[jj][i];
          v.s!(jj,j,x*c+z*s); v.s!(jj,i,z*c-x*s)
        end # for jj in 1..n
        z = f.hypot(h)
        w[j] = z
        if z.nonzero?
          z = one/z; c = f * z; s = h * z
        end # if z.nonzero?
        f=c*g+s*y; x=c*y-s*g;
        for jj in 1..m
          y=a[jj][j]; z=a[jj][i]
          a.s!(jj,j,y*c+z*s); a.s!(jj,i,z*c-y*s)
        end # for jj in 1..m
      end # for j in l..nm
      rv1[l] = zero
      rv1[k] = f
      w[k] = x
    end # for its in 1..iter_max
    k -= 1
  end # while k >=1
end

.test(n0, verbose = true, complex = false) ⇒ Object

Consistency test of class Mat.

# File 'lib/linalg.rb', line 1105

def Mat.test(n0, verbose = true, complex = false )
  puts "Doing Mat.test( n = #{n0}, verbose = #{verbose}," +
  " complex = #{complex} ) for R.prec = #{R.prec}:"
  puts "******************************************************"
  t1 = Time.now
  s = R.c0
  puts "class of s is " + s.class.to_s
  i = n0 + 137
  a = Mat.tob(n0, i, complex) 
  i += 1
  b = Mat.tob(n0, i, complex) 
  i += 1
  c = Mat.tob(n0, i, complex) 
  x = Vec.tob(n0, i, complex) 
  i += 1
  y = Vec.tob(n0, i, complex) 
  i += 1
  s1 = complex ? C.ran(i) : R.ran(i)
  i += 1
  s2 = complex ? C.ran(i) : R.ran(i)

  unit = a.to_1

  a0 = a.clone
  b0 = b.clone
  c0 = c.clone

  abc0 = a0.abs + b0.abs + c0.abs

  ac = a.clone
  a = b
  r = a
  l = b
  ds = r.dis(l)
  puts "assignment of variables: ds = " + ds.to_s if verbose
  s += ds

  a = ac
  r = a
  l = ac
  ds = r.dis(l)
  puts "assignment of variables 2: ds = " + ds.to_s if verbose
  s += ds
 
  r = (a + b) + c
  l = a + (b + c)
  ds = r.dis(l)
  puts "associativity of +: ds = " + ds.to_s if verbose
  s += ds

  r = (a - b) + c
  l = a - (b - c)
  ds = r.dis(l)
  puts "associativity of -: ds = " + ds.to_s if verbose
  s += ds

  r = (a * b) * c
  l = a * (b * c)
  ds = r.dis(l)
  puts "associativity of *: ds = " + ds.to_s if verbose
  s += ds

  r = (a + b) * s1
  l = a * s1 + b * s1
  ds = r.dis(l)
  puts "distributivity of multiplication by scalars: ds = " + 
    ds.to_s if verbose
  s += ds

  r = c * (s1*s2)
  l = (c * s1) * s2
  ds = r.dis(l)
  puts "distributivity of multiplication by scalars: ds = " + 
    ds.to_s if verbose
  s += ds

  r = a
  l = -(-a)
  ds = r.dis(l)
  puts "idempotency of unary minus: ds = " + ds.to_s if verbose
  s += ds

  r = (a + b).spr(c)
  l = a.spr(c) + b.spr(c)
  ds = r.dis(l)
  puts "distributivity of spr: ds = " + ds.to_s if verbose
  s += ds

  r = (a + b) * c
  l = a * c + b * c
  ds = r.dis(l)
  puts "distributivity of matrix multiplication: ds = " + 
    ds.to_s if verbose
  s += ds

  r = (a * b) * x
  l = a * (b * x)
  ds = r.dis(l)
  puts "action on vectors 1: ds = " + ds.to_s if verbose
  s += ds

  r = (a + b) * x
  l = a * x + b * x
  ds = r.dis(l)
  puts "action on vectors 2: ds = " + ds.to_s if verbose
  s += ds

  r = b * (x + y)
  l = b * x + b * y
  ds = r.dis(l)
  puts "action on vectors 3: ds = " + ds.to_s if verbose
  s += ds

  r = c * (x * s1)
  l = (c * s1) * x 
  ds = r.dis(l)
  puts "action on vectors 4: ds = " + ds.to_s if verbose
  s += ds

  if complex == false
    r = unit
    l = a * a.pseudo_inv
    ds = r.dis(l)
    puts "pseudo inverse is right inverse: ds = " + ds.to_s if verbose
    s += ds

    r = unit
    l = a.pseudo_inv * a
    ds = r.dis(l)
    puts "pseudo inverse is left inverse: ds = " + ds.to_s if verbose
    s += ds
  else
    puts "test of pseudo inverse left out, since not implemented for complex"
  end

  aMem = a.clone
  bMem = b.clone
  cMem = c.clone

# Testing the access functions, under harsh conditions with
# inserting double transposition

  for i in 1..a.dim1
    for j in 1..a.dim2
      b.s!(i,j,a[i][j]) # copies a to b
    end
  end
  
  l = a
  r = a.trp.trp
  ds = r.dis(l)
  puts "test of double transposition: ds = " + ds.to_s if verbose
  s += ds

  for i in 1..b.dim1
    for j in 1..b.dim2
      c.s!(i,j,b[i][j]) # copies b to c
    end
  end
# Finally c should have the content of a

  r = a
  l = c
  ds = r.dis(l)
  puts "test of access functions: ds = " + ds.to_s if verbose
  s += ds

  a = aMem; b = bMem; c = cMem

  abc = a.abs + b.abs + c.abs

  l = abc
  r = abc0
  ds = r.dis(l)
  puts "heavy test of assignment: ds = " + ds.to_s if verbose

  t2 = Time.now

  if verbose
    puts
    a.prn("a")
    puts
    b.prn("b")
    puts
    c.prn("c")
    puts
    s1.prn("s1")
    puts
    s2.prn("s2")
    puts
    x.prn("x")
    puts
    y.prn("y")
  end

  puts "class of s is " + s.class.to_s + " ."
  puts "The error sum s is " + s.to_s + " ."
  puts "It should be close to 0."
  puts "Computation time was " + (t2-t1).to_s
  s
end

.tob(n, i, complex = false) ⇒ Object

Generates a test object, here a n times n matrix with random elements. This object depends rather chaotically on the integer parameter i. If the last argument is ‘false’ the test matrix will have R-typed elements, and C-typed elements else.

# File 'lib/linalg.rb', line 506

def Mat.tob(n,i, complex = false)
  if complex
    ri = C.tob(i)
    zero = ri.to_0
    res=Mat.new(n, n, zero)
    rg1 = Ran.new(-ri.re,ri.re)
    rg2 = Ran.new(-ri.im,ri.im)
    for j in 1..n
      for k in 1..n
        yjk = C.new(rg1.ran,rg2.ran)
        res.s!(j,k,yjk)
      end
    end
    res
  else
    ri = R.tob(i)
    zero = ri.to_0
    res=Mat.new(n, n, zero)
    rg = Ran.new(-ri,ri)
    for j in 1..n
      for k in 1..n
        yjk = rg.ran
        res.s!(j,k,yjk)
      end
    end
    res
  end
end

Instance Method Details

#*(v) ⇒ Object

Multiplication of a Mat with either a Mat, Vec, or Numeric

# File 'lib/linalg.rb', line 926

def *(v)
  d1 = dim1
  fail "dim1 == 0" if d1.zero?
  d2 = dim2
  fail "dim2 == 0" if d2.zero?
  zero=@x[0].x[0].to_0
  if v.is_a?(Mat)
    d3 = v.dim1
    fail "dimenson mismatch" unless d3 == d2
    d4 = v.dim2
    fail "dim4 == 0" if d4.zero?
  # we better produce the d1 row-vectors in turn
    a = Array.new(d1)
    for i in 0...d1
      vi = Vec.new(d4,zero)
      for j in 0...d4
        vij = zero
        for k in 0...d2
          vij += @x[i].x[k] * v.x[k].x[j]
        end
        vi.x[j] = vij
      end
      a[i] = vi
    end
    res = Mat.new(a)
  elsif v.is_a?(Vec) 
    d3 = v.dim
    fail "dimenson mismatch" unless d3 == d2
    res = Vec.new(d1,zero)
    for i in 0...d1
      vi = zero
      for j in 0...d2
        vi += @x[i].x[j] * v.x[j]
      end
      res.x[i] = vi
    end
  elsif v.is_a?(Numeric) # multiplication with scalar
    res = clone
    for i in 1..d1
      res[i] *= v
    end
  else
    fail "can't multiply with this object"
  end
  res
end

#+(v) ⇒ Object

Returns self + v, where v is a Mat

# File 'lib/linalg.rb', line 793

def +(v)
  fail "Object can't be added to a Mat." unless v.is_a?(Mat)
  fail "Dimension mismatch." unless dim == v.dim
  res = clone
  for i in 1..dim
    res[i] += v[i]
  end
  res
end

#-(v) ⇒ Object

Returns self - v , where v is a Mat

# File 'lib/linalg.rb', line 804

def -(v)
  fail "Object can't be subtracted from a Mat." unless v.is_a?(Mat)
  fail "Dimension mismatch." unless dim == v.dim
  res = clone
  for i in 1..dim
    res[i] -= v[i]
  end
  res
end

#-@ ⇒ Object

Unary minus operator. Returns - self.

# File 'lib/linalg.rb', line 784

def -@
  res = clone
  for i in 1..dim
    res[i] = -res[i]
  end
  res
end

#<=>(v) ⇒ Object

The order relation is here lexicographic ordering of lists. Needed only for book-keeping purposes.

# File 'lib/linalg.rb', line 975

def <=> (v)
  d1 = dim; d2 = v.dim
  if d1 < d2
    return -1
  elsif d1 > d2 
    return 1
  else
    for i in 1..d1
      ci = self[i] <=> v[i]
      return ci unless ci == 0
    end
  end
  return 0
end

#[](i) ⇒ Object

Reading row vectors. Valid indexes are 1,…,dim1.

# File 'lib/linalg.rb', line 741

def [](i)
  @x[i-1]
end

#[]=(i, a) ⇒ Object

Setting row vectors. Valid indexes are 1,…,dim1. Notice that setting matrix elements via [][]= is not permanent. For setting matrix elements the method s! is provided.

# File 'lib/linalg.rb', line 748

def []=(i,a) 
  @x[i-1] = a
end

#abs ⇒ Object

Absolute value, always real

# File 'lib/linalg.rb', line 1011

def abs
  if complex?
    abs2.re.sqrt
  else
    abs2.sqrt
  end
end

#abs2 ⇒ Object

Square of absolute value.

# File 'lib/linalg.rb', line 1006

def abs2
  spr(self)
end

#clone ⇒ Object

Returns an independent copy of self.

# File 'lib/linalg.rb', line 497

def clone
  Mat.new(self)
end

#complex? ⇒ Boolean

Returns:

• (Boolean)

# File 'lib/linalg.rb', line 1034

def complex? 
  return nil if dim.zero?
  @x[0].complex?
end

#conj ⇒ Object

Returns the Hermitian conjugate of self.

# File 'lib/linalg.rb', line 835

def conj
  d1 = dim1
  d2 = dim2
  fail "dim1 == 0" if d1.zero?
  fail "dim2 == 0" if d2.zero?
  zero = @x[0].x[0].to_0
  res = Mat.new(d2,d1,zero)
  for i in 0...d1
    for j in 0...d2
      sij = @x[i].x[j]
      res.s!(j+1,i+1,sij.conj)
    end
  end
  res
end

#dim ⇒ Object

Returns the ‘dimension’ of the matrix, i.e. the number of its row-vectors. This thus is m for a ‘m times n - matrix’.

# File 'lib/linalg.rb', line 400

def dim; @x.size; end

#dim1 ⇒ Object

Let self be a (m,n)-matrix (also called a m times n matrix) then dim1 == m

# File 'lib/linalg.rb', line 404

def dim1; @x.size; end

#dim2 ⇒ Object

Let self be a (m,n)-matrix (also called a m times n matrix) then dim2 == n

# File 'lib/linalg.rb', line 408

def dim2
  return 0 if dim.zero?
  self[1].dim
end

#dis(v) ⇒ Object

Relative distance between matrices

# File 'lib/linalg.rb', line 1020

def dis(v)
  a = abs
  b = v.abs
  d = (self - v).abs
  s = a + b
  return R.c0 if s.zero?
  d1 = d/s
  d < d1 ? d : d1
end

#each ⇒ Object

# File 'lib/linalg.rb', line 990

def each
  @x.each{ |c| yield c}
end

#inv ⇒ Object

In most cases (‘up to a subst of measure zero’) the pseudo-inverse is also the inverse.

# File 'lib/linalg.rb', line 1100

def inv
  pseudo_inv
end

#prn(name) ⇒ Object

Prints the content of self and naming the output.

# File 'lib/linalg.rb', line 775

def prn(name)
  for i in 1..dim1 
    for j in 1..dim2
      puts " #{name}[#{i}][#{j}] = " + self[i][j].to_s
    end
  end
end

#pseudo_inv(acc = 0) ⇒ Object

Returns the pseudo-inverse (als known as Penrose inverse) of self. If the argument acc is not zero, the discontinous treatment of singular values near zero is replaced by a continuous one. Notice that the pseudo inverse always exists, and that the pseudo- inverse of a (m,n)-matrix is a (n,m)-matrix. If the argument acc is not zero, the pseudo-inverse is the first component of a 4-array res, which also contains the the intermediary quantities a, w, v resulting from the call

Mat.svdcmp(a,w,v). a == res[1], w == res[2], v == res[3].

Especially the list w of the original singular values is thus made accessible, so that one can judge whether their processing controlled by the parameter acc was reasonable. The pseudo_inverse is the most useful and stable mehod to solve linear equations: Let a be a (m,n)-matrix, and b a m-vector. The equation

a * x = b (i)

determines a n-vector x as

x = a.pseudo_inv * b

which is a solution of (i) if there is one. If there are many solutions it is the one of minimum absolute value, and if there is no solution it comes closest to be a solution: It minimizes the defect

(a * x - b).abs

Simply great! No need for LU-decompositions any more.

# File 'lib/linalg.rb', line 1063

def pseudo_inv(acc = 0)
  if complex?
    fail "Pseudo inverse not yet impemented for complex matrices"
  end
  m = dim1
  fail "dim1 == 0" if m.zero?
  n=dim2
  fail "dim2 == 0" if n.zero?
  rightframe = m >= n
  a = clone
  a = a.trp if rightframe == false
  m = a.dim1;  n = a.dim2
  zero = a[1][1].to_0
  v = Mat.new(n,n,zero)
  w = Vec.new(n,zero)
  vr = Mat.new(n,m,zero) #  sic
  Mat.svdcmp(a,w,v)
  wi = w.pseudo_inv(acc)  # w does not come out orderd
  for i in 1..n
    for j in 1..m
      sum = zero
      for k in 1..n
        sum += v[i][k] * wi[k] * a[j][k]
      end
      vr.s!(i,j,sum)
    end
  end
  vr = vr.trp if rightframe == false
  if acc.zero? 
    vr
  else
    [vr,a,w,v]
  end
end

#s!(i, j, a) ⇒ Object

The ‘s’ stands for ‘set (value)’ and the ‘!’ this is a method by which self changes (non-constant or mutating method).

# File 'lib/linalg.rb', line 755

def s!(i,j,a)
  si = Vec.new(self[i]) # don't change the row-vector self[i]
    # itself. Such changes are subject to subtle side effects.
    # Worcing on a copy is safe.
  si[j] = a   # changing si here is normal syntax 
  self[i] = si # this is OK, of course 
    # self[i] = Vec.new(si) would work too, but would cause more work 
end

#spr(v) ⇒ Object

Scalar product of matrices.

# File 'lib/linalg.rb', line 995

def spr(v)
  fail "dimension mismatch" unless dim == v.dim
  return nil if dim.zero?
  s = self[1].spr(v[1]) 
  for i in 2..dim
    s += self[i].spr(v[i]) 
  end
  s
end

#square? ⇒ Boolean

Returns:

• (Boolean)

# File 'lib/linalg.rb', line 1030

def square?
  dim1 == dim2
end

#to_0 ⇒ Object

‘to zero’ Returns a matrix with he same dimensions as self, but with all matrix elements set to zero.

# File 'lib/linalg.rb', line 855

def to_0
  d1 = dim1
  d2 = dim2
  fail "dim1 == 0" if d1.zero?
  fail "dim2 == 0" if d2.zero?
  zero = @x[0].x[0].to_0
  Mat.new(d2,d1,zero)
end

#to_1 ⇒ Object

‘to one’. Returns a matrix with he same dimensions as self, but with all matrix elements set to 1.

# File 'lib/linalg.rb', line 866

def to_1
  d1 = dim1
  d2 = dim2
  fail "dim1 == 0" if d1.zero?
  fail "dim2 == 0" if d2.zero?
  fail "dim1 != dim2" unless d1 == d2
  unit = @x[0].x[0].to_1
  diag = Vec.new(d1,unit)
  Mat.new(diag)
end

#to_s ⇒ Object

Returns a string which consists of a list of the strings which represent the row-vectors.

# File 'lib/linalg.rb', line 766

def to_s
  res = "\n Mat"
  for i in 0...dim
    res += "\n    " + x[i].to_s
  end
  res + "\n end Mat"
end

#trp ⇒ Object

Returns the transposed of self.

# File 'lib/linalg.rb', line 815

def trp # implementation without s!
  d1 = dim1
  d2 = dim2
  fail "dim1 == 0" if d1.zero?
  fail "dim2 == 0" if d2.zero?
  zero = @x[0][0].to_0
  # self has d1 row-vectors of length d2. 
  # The result of transposing has d2 row-vectors of length d1.
  v = Array.new(d2)
  for j in 0...d2
    vj = Vec.new(d1,zero)
    for i in 0...d1
      vj.x[i] = @x[i].x[j]
    end
    v[j] = vj
  end
  Mat.new(v)
end