Class: NMatrix

Inherits:
Object
  • Object
show all
Defined in:
lib/mixed_models/nmatrix_methods.rb

Overview

Copyright © 2015 Alexej Gossmann

Class Method Summary collapse

Instance Method Summary collapse

Class Method Details

.block_diagonal(*params) ⇒ Object

Generate a block-diagonal NMatrix from the supplied 2D square matrices.

Arguments

  • params - An array that collects all arguments passed to the method. The method can receive any number of arguments. The last entry of params is a hash of options from NMatrix#initialize. All other entries of params are the blocks of the desired block-diagonal matrix, which are supplied as square 2D NMatrix objects.

Usage

a = NMatrix.new([2,2],[1,2,3,4])
b = NMatrix.new([1,1],[123],dtype: :int32)
c = NMatrix.new([3,3],[1,2,3,1,2,3,1,2,3], dtype: :float64)
m = NMatrix.block_diagonal(a,b,c,dtype: :int64, stype: :yale)
    => 
    [
      [1, 2,   0, 0, 0, 0]
      [3, 4,   0, 0, 0, 0]
      [0, 0, 123, 0, 0, 0]
      [0, 0,   0, 1, 2, 3]
      [0, 0,   0, 1, 2, 3]
      [0, 0,   0, 1, 2, 3]
    ]



134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
 
# File 'lib/mixed_models/nmatrix_methods.rb', line 134

def block_diagonal(*params)
  options = params.last.is_a?(Hash) ? params.pop : {}

  block_sizes = [] #holds the size of each matrix block
  params.each do |b|
    raise ArgumentError, "Only NMatrix objects allowed" unless b.is_a?(NMatrix)
    raise ArgumentError, "Only 2D matrices allowed" unless b.shape.size == 2
    raise ArgumentError, "Only square matrices allowed" unless b.shape[0] == b.shape[1]
    block_sizes << b.shape[0]
  end

  bdiag = NMatrix.zeros(block_sizes.sum, options)
  (0...params.length).each do |n|
    # First determine the size and position of the n'th block in the block-diagonal matrix
    k = block_sizes[n]
    block_pos = block_sizes[0...n].sum
    # populate the n'th block in the block-diagonal matrix
    (0...k).each do |i|
      (0...k).each do |j|
        bdiag[block_pos+i,block_pos+j] = params[n][i,j]
      end
    end
  end

  return bdiag
end

Instance Method Details

#khatri_rao_rows(mat) ⇒ Object

Compute a simplified version of the Khatri-Rao product of self and other NMatrix mat. The i’th row of the resulting matrix is the Kronecker product of the i’th row of self and the i’th row of mat.

Arguments

  • mat - A 2D NMatrix object

Usage

a = NMatrix.new([3,2], [1,2,1,2,1,2], dtype: dtype, stype: stype)
b = NMatrix.new([3,2], (1..6).to_a, dtype: dtype, stype: stype)
m = a.khatri_rao_rows b # =>  [ [1.0, 2.0,  2.0,  4.0]
                                [3.0, 4.0,  6.0,  8.0]
                                [5.0, 6.0, 10.0, 12.0] ]

Raises:

  • (NotImplementedError) 


43
44
45
46
47
48
49
50
51
52
53
54
55
 
# File 'lib/mixed_models/nmatrix_methods.rb', line 43

def khatri_rao_rows(mat)
  raise NotImplementedError, "Implemented for 2D matrices only" unless self.dimensions==2 and mat.dimensions==2
  n = self.shape[0]
  raise NotImplementedError, "Both matrices must have the same number of rows" unless n==mat.shape[0]
  m = self.shape[1]*mat.shape[1]
  prod_dtype = NMatrix.upcast(self.dtype, mat.dtype)
  khrao_prod = NMatrix.new([n,m], dtype: prod_dtype)
  (0...n).each do |i|
    kronecker_prod = self.row(i).kron_prod_1D mat.row(i)
    khrao_prod[i,0...m] = kronecker_prod
  end
  return khrao_prod
end

#kron_prod_1D(v) ⇒ Object

Compute the the Kronecker product of two row vectors (NMatrix of shape [1,n])

Arguments

  • v - A NMatrix of shape [1,n] (i.e. a row vector)

Usage 

a = NMatrix.new([1,3], [0,1,0])
b = NMatrix.new([1,2], [3,2])
a.kron_prod_1D b #  =>  [ [0, 0, 3, 2, 0, 0] ]



17
18
19
20
21
22
23
24
25
 
# File 'lib/mixed_models/nmatrix_methods.rb', line 17

def kron_prod_1D(v)
  unless self.dimensions==2 && v.dimensions==2 && self.shape[0]==1 && v.shape[0]==1
    raise ArgumentError, "Implemented for NMatrix of shape [1,n] (i.e. one row) only."
  end
  #TODO: maybe some outer product function from LAPACK would be more efficient to compute for m
  m = self.transpose.dot v
  l = self.shape[1]*v.shape[1]
  return m.reshape([1,l])
end

#matrix_valued_solve(rhs_mat) ⇒ Object

Solve a matrix-valued linear system A * X = B, where A is a nxn matrix, B and X are nxp matrices.

Arguments

  • rhs_mat - the 2D matrix on the right hand side of the linear system

Usage

a = NMatrix.random([10,10], dtype: :float64)
b = NMatrix.random([10,5], dtype: :float64)
x = a.solve(b)



69
70
71
72
73
74
75
76
77
 
# File 'lib/mixed_models/nmatrix_methods.rb', line 69

def matrix_valued_solve(rhs_mat)
  # This is a workaround. See NMatrix issue:
  # https://github.com/SciRuby/nmatrix/issues/374
  rhs_mat_t = rhs_mat.transpose
  lhs_mat = self.clone
  NMatrix::LAPACK::clapack_gesv(:row, lhs_mat.shape[0], rhs_mat.shape[1], 
                                lhs_mat, lhs_mat.shape[0], rhs_mat_t, rhs_mat.shape[0])
  return rhs_mat_t.transpose
end

#triangular_solve(uplo, rhs) ⇒ Object

Solve a linear system A * X = B, where A is a lower triangular matrix, X and B are vectors or matrices.

Arguments

  • uplo - flag indicating whether the matrix is lower or upper triangular; possible values are :lower and :upper

  • rhs - the right hand side, an NMatrix object

Usage

a = NMatrix.new(3, [4, 0, 0, -2, 2, 0, -4, -2, -0.5], dtype: :float64)
b = NMatrix.new([3,1], [-1, 17, -9], dtype: :float64)
x = a.triangular_solve(:lower, b)
a.dot x # => [ [-1.0]   [17.0]   [-9.0] ]

Raises:

  • (ArgumentError) 


95
96
97
98
99
100
101
102
103
104
 
# File 'lib/mixed_models/nmatrix_methods.rb', line 95

def triangular_solve(uplo, rhs)
  raise(ArgumentError, "uplo should be :lower or :upper") unless uplo == :lower or uplo == :upper
  b = rhs.clone
  # this is the correct function call; it came up in during
  # discussion in https://github.com/SciRuby/nmatrix/issues/374
  NMatrix::BLAS::cblas_trsm(:row, :left, uplo, false, :nounit, 
                            b.shape[0], b.shape[1], 1.0, self, self.shape[0],
                            b, b.shape[1])
  return b
end