Monday, July 22, 2019

Multivariate Normal Sampling with Scala and ND4j

Target audience: Intermediate
Estimated reading time: 6'

The multivariate normal sampling function plays a pivotal role in several machine learning methods, including Markov chains Monte Carlo (MCMC) and the contextual multi-armed bandit. The foundation of this sampling is data vectorization, a technique familiar to both Python developers and data scientists.
In this article, we outline the process of implementing multivariate normal sampling with the aid of ND4j.

Table of contents


 Notes: 
  • The post requires some knowledge of data vectorization (numpy, datavec, ND4j..) as well as Scala programming language.
  • The code associated with this article is written using Scala 2.11.8

Vectorization

Python Numpy is a well-known and reliable vectorized linear algebra library which is a foundation of scientific (SciPy) and machine learning (Sciktlearn) libraries. No serious data scientific projects can reasonably succeeds with the power and efficiency of numpy. 
The vectorization of datasets is the main reason behind the performance of machine learning models (training and prediction) build in Python.

 Is there a similar linear algebra library, supporting vectorization, available to Scala and Spark developers? Yes, ND4j

ND4j, BLAS and LAPACK

ND4j library replicates the functionality of numpy for Java developers. ND4j can be downloaded as an independent library or as component of the deep learning library, deeplearning4j. It leverages the BLAS and LAPACK libraries.
From a Java developer perspective, the data represented by an NDArray is stored outside of the Java Virtual Machine. This design has the following benefits:
  • Avoid taxing the Garbage Collector
  • Interoperability with high-performance BLAS libraries such as OpenBlas
  • Allow number of array elements exceeds Int.MaxValue
The BLAS (Basic Linear Algebra Subprograms) are functions performing basic vector and matrix operations. The library is divided in 3 levels:
  • Level 1 BLAS perform scalar, vector and vector-vector operations,
  • Level 2 BLAS perform matrix-vector operations
  • Level 3 BLAS perform matrix-matrix operations. 
OpenBLAS is an optimized Basic Linear Algebra Subprograms (BLAS) library based on GotoBLAS2,  a C-library of linear algebra supporting a large variety of micro-processor. Its usage is governed by the BSD license.
LAPACK are Fortran routines for solving systems of simultaneous linear equations, least-squares solutions of linear systems of equations, eigenvalue problems, and singular value problems and matrix factorizations.

Implicit ND4j array conversion

The first step is to create a implicit conversation between ND4j and Scala data types.  
The following code convert an array of double into a INDArray using org.nd4j.linalg.factory.Nd4j java class and its constructor create(double[] values, int[] shape)
  • In case of a vector, the shape is defined in Scala as  Array[Int](size_vector)
  • In case of a matrix, the shape is defined as Array[Int](numRows, numCols). 
The following snippet implement a very simple conversion from a Scala array to a INDArray



@throws(clazz = classOf[IllegalArgumentException])
 
implicit def double2NDArray(values: Array[Double]): INDArray = {
  require(values.nonEmpty, "ERROR: ND4, conversion ...")

  Nd4j.create(values, Array[Int](1, values.length))
}











Multivariate Normal distribution

The sampling of a multivariate normal distribution is defined by the formula 

\[\mathbb{N}(\mu, \Sigma )=\frac{1}{\sqrt{(2\pi)^n \left | \Sigma \right |}} e^{\frac{1}{2}(x-\mu)^T {\Sigma '}^{-1}(x-\mu)}\] A generalized definition adds a very small random perturbation factor r |r| <= 1 on the variance value (diagonal of the covariance matrix)
\[\Sigma ' = \Sigma + r.I\] Sigma is the covariance matrix and the mu is the mean value. 

 The computation of the multivariate normal sampling can be approximated by the Cholesky decomposition. In a nutshell, the Cholesky algorithm decompose a positive-definite matrix into a product of two matrix


    

  • lower triangle matrix
    • 

  • transposition of its conjugate
    • 


It serves the same purpose as the ubiquitous LU decomposition with less computation. \[\mathbb{N}(\mu, \Sigma) \sim \mu + Chol(\Sigma).Z^T\] where \[L=Chol(\Sigma)\] and \[L.L^T=\Sigma\]. The vector Z is a multivariate normal noise \[Z= \{ z_i | z_i=N(0, 1)\}\]

The following implementation relies on the direct invocation of LAPACK library function potrf. The LAPACK functionality are accessed through the BLAS wrapper getBlasWrapper.












 final def choleskyDecomposition(matrix: INDArray): INDArray = {
     val matrixDup = matrix.dup
     Nd4j.getBlasWrapper.lapack.potrf(matrixDup, false)
     matrixDup
 }







Note that the matrix is duplicated prior to the LAPACK function call as we do not want to alter the input matrix. 

Let's implement the multivariate Normal sampler with perturbation using the Cholesky decomposition.







 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45


@throws(clazz = classOf[IllegalArgumentException])
@throws(clazz = classOf[IllegalStateException])
final def multiVariateNormalSampling(
   mean: INDArray,
   cov: INDArray,
   perturbation: Double = 0.0): INDArray = {

    import scala.util.Random
    require(cov.size(0) == mean.length, s"Sigma shape ${cov.size(0)} should be mean size ${mean.length}")
    require(cov.size(1) == mean.length, s"Sigma shape ${cov.size(1)} should be ${mean.length}")

    try {
      // Add a perturbation epsilon value to the covariance matrix
      // if it is defined
      val perturbMatrix =
         if(perturbation > 0.0)
            cov.add( squareIdentityMatrix(cov.size(0), 
                           perturbation*(Random.nextDouble-0.5)))
         else
            cov

      // Apply the Cholesky decomposition
      val perturbed: INDArray = choleskyDecomposition(perturbMatrix)
      // Instantiate a normal distribution
      val normalPdf = new NormalDistribution(
           new DefaultRandom, 
           0.0, 
           1.0)

      val normalDensity = Array.fill(mean.size(0))(normalPdf.sample)
      val normalNDArray: INDArray = normalDensity

       // This is an implementation of the Dot product
      val normalCov = perturbed.mmul(normalNDArray.transpose)
      // Add the normalized perturbed covariance to the mean value
      mean.add(normalCov)
   }
   catch {
       case e: org.nd4j.linalg.api.blas.BlasException =>
           throw new IllegalStateException(e.getMessage)
       case e: org.apache.commons.math3.exception.NotStrictlyPositiveException =>
           throw new IllegalStateException(e.getMessage)
       case e: Throwable =>
           throw new IllegalStateException(e.getMessage)
   }
}









Let's look at the full implementation of the multi-variate normal sampling.

The first step validates the shape of the mean and covariance input matrices [line 8, 9]. As mentioned earlier, the generalized normal distribution introduces an optional random perturbation of small magnitude (~1e-4) [line 14-17] that is useful for application that requires some stochastic

The 'perturbed' covariance matrix is factorized using the Cholesky decomposition [line 22]. The normal probability density function (mean 0.0 and standard deviation 1.0) is used to generate random values [line 24-30] which is then applied to the covariance matrix [line 33].

The normal randomized variance is added to the vector of mean values [line 35].



For the sake of convenience, the multivariate normal density function uses the Apache Math common 3.5 API [line 24].





References




---------------------------
Patrick Nicolas has over 25 years of experience in software and data engineering, architecture design and end-to-end deployment and support with extensive knowledge in machine learning. 
He has been director of data engineering at Aideo Technologies since 2017 and he is the author of "Scala for Machine Learning" Packt Publishing ISBN 978-1-78712-238-3












Posted by






No comments:
  






















Labels:
,
,
,
,
,
,
















          

        

          

        

Sunday, May 12, 2019


          

        










Genetic Algorithms III: Solver











Target audience: Advanced
Estimated reading time: 9'

Genetic algorithms fall under the umbrella of optimization algorithms, which are specifically designed to identify the peak or trough of a given function.
In this article, we present the concluding segment of our series on the Scala-based implementation of genetic algorithms, focusing on the reproduction cycle and the solver.


Table of contents


Follow me on LinkedIn


Introduction


In the first post on genetic algorithms, you learned about the key elements of genetic algorithms (Chromosomes, Genes and Population). Genetic Algorithms I: Foundation


This second part introduces the concept and implementation of genetic operations (cross-over, mutation and selection) on a population of chromosomes. Genetic Algorithms II: Operators
This 3rd and final post on genetic algorithms, explores the application of genetic algorithm as a solver or optimizer




Note: 
For the sake of readability of the implementation of algorithms, all non-essential code such as error checking, comments, exception, validation of class and method arguments, scoping qualifiers or import is omitted.





Reproduction cycle




Let's wrap the reproduction cycle into a Reproduction class that uses the scoring
function score. The class defines a random generator (line 2) used in identifying the cross-over and mutation rates. The key method in the reproduction cycle is mate (line 4).






1
2
3
4
5
6
7
8
class Reproduction[T <: Gene](score: Chromosome[T] => Unit) {     
  val rand = new Random(System.currentTimeMillis)

  def mate(
      population: Population[T], 
      config: GAConfig, 
      cycle: Int): Boolean = {}
}









The reproduction function, mate, implements the sequence of the three
genetic operators as a workflow:

select for the selection of chromosomes from the population (line 8)

+- for the crossover between chromosomes (line 9)

^ for the mutation of the chromosomes (line 10).







 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
def mate(
    population: Population[T], 
    config: GAConfig, 
    cycle: Int): Boolean = population.size match {

  case 0 | 1 | 2 => false
  case _ => {
      population.select(score, config.softLimit(cycle))
      population +- (1.0 - Random.nextDouble*config.xover)
      population ^ (1.0 - Random.nextDouble*config.mu)
      true
  }
}










This method returns true (line 11) if the size of the population is larger than 2 (line 6). The last
element of the puzzle (reproduction cycle) is the exit condition. There are two options for estimating
that the reproducing cycle is converging:


    

  • Greedy: In this approach, the objective is to evaluate whether the n fittest
chromosomes have not changed in the last m reproduction cycles
    • 

  • Loss (or objective) function: This approach is similar to the convergence criteria for the
training of supervised learning.


Solver




The last class GASolver manages the reproduction cycle and evaluates the exit
condition or the convergence criteria:







 1
 2
 3
 4
 5
 6
 7
 8
 9
10
class GASolver[T <: Gene](
    config: GAConfig, 
    score: Chromosome[T] => Unit
) extends PipeOperator[Population[T], Population[T]] {
  
    var state: GAState = GA_NOT_RUNNING
 
    def solve(in: Population[T]): Population[T]) {}
     ...
}









The class GASolver implements the data transformation, solve which transforms a population to
another one, given a configuration of the genetic algorithm, config (line 2) and a scoring method, score (line 3). The class is also responsible for maintaining and updating the state of the reproduction cycle (line 6).






 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
def solve(pop: Population[T]): Population[T] = 
    if(in.size > 1) {
       val reproduction = Reproduction[T](score)
       state = GA_RUNNING
     
       Range(0, config.maxCycles).find(n => { 
         reproduction.mate(pop, config, n) match {
             case true => converge(pop, n) != GA_RUNNING
             case false => {}
         }
       }) match {
           case Some(n) => pop
           case None => {}
       }
    }
    else
       pop





As mentioned previously, the solve method transforms a population by applying the sequence of genetic operators.


It instantiates a new reproduction cycle (line 3) and set the internal state of the genetic algorithm as RUNING (line 4). The mating of all the chromosomes in the population is implemented iteratively (lines 6 - 10). It exits when either the maximum number of cycles (line 6) is reached or the reproduction converged to a single solution (line 8).




Putting all together
Let's apply the genetic solver to selecting a strategy (i.e. Trading, Marketing,...). The steps in the execution of the solver (or optimizer) using the genetic algorithm are:


    

  1. Initialize the execution & genetic algorithm configuration parameters (lines 1-6).
    2. 

  2. Define a soft limit on the maximum size of the chromosomes population (line 7)
    3. 

  3. Specify the quantization or discretization function (line 9)
    4. 

  4. Define the scoring function applied to a chromosome if a trading signal as gene (lines 15-17)
    5. 

  5. Initialize the population of chromosomes as trading strategies (line 20)
  6. Initialize the solver (lines 23 & 24)
    7. 

  7. execute the genetic algorithm as a iterative sequence of reproduction cycles (line 28)





 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
val XOVER = 0.8  // Probability or ratio for cross-over
val MU = 0.4         // Probability or ratio for mutation
val MAX_CYCLES = 400  // Maximum number of iterations during the optimization
val CUTOFF_SLOPE = -0.003  // Slope for the linear soft limit
val CUTOFF_INTERCEPT = 1.003 // Intercept value for the linear soft limit
val R = 1024  // quantization ratio for conversion Int <-> Double
val softLimit = (n: Int) => CUTOFF_SLOPE*n + CUTOFF_INTERCEPT    

implicit val digitize = new Discretization(R)


  // Define the scoring function for the chromosomes (i.e. Trading 
  // strategies) as the sum of the score of the genes 
  // (i.e. trading signals) in this chromosome (i.e strategy)
val scoring = (chr: Chromosome[Signal]) =>  {
    val signals: List[Gene] = chr.code
    chr.unfitness = signals.foldLeft(0.0)((sum, s) => sum + s.score)
}

val population = Population[Signal]((strategies.size <<4), strategies)

// Configure, instantiates the GA solver for trading signals
val config = GAConfig(XOVER, MU, MAX_CYCLES, softLimit)
val gaSolver = GASolver[Signal](config, scoring)

  // Extract the best population and the fittest chromosomes = 
  // trading strategies from this final population.
val bestPopulation = gaSolver solve population

bestPopulation.fittest(1)
     .getOrElse(ArrayBuffer.empty)
     .foreach(
          ch => logger.info(s"$name Best strategy: ${ch.symbolic("->")}")
      )









 Thank you for reading this article. For more information ...


References


---------------------------
Patrick Nicolas has over 25 years of experience in software and data engineering, architecture design and end-to-end deployment and support with extensive knowledge in machine learning. 
He has been director of data engineering at Aideo Technologies since 2017 and he is the author of "Scala for Machine Learning" Packt Publishing ISBN 978-1-78712-238-3










Posted by






No comments:
  






















Labels:
,
,
,
,
,
















          

        

          

        

Wednesday, April 3, 2019


          

        










Genetic Algorithms II: Operators











Target audience: Advanced
Estimated reading time: 6'

In the realms of science and engineering, genetic algorithms serve dual purposes: as adaptable algorithms addressing real-world challenges and as computational representations of natural evolutionary mechanisms. 

This article stands as the subsequent chapter in our series on genetic algorithms. Within, we detail the Scala-based implementation of genetic operators, encompassing selection, cross-over, and mutation, acting upon a chromosome population.

Table of contents


Follow me on LinkedIn
 
Introduction







In the first article on genetic algorithms, you learned about the key elements of genetic algorithms.Genetic Algorithms I: Foundation


    

  • Chromosomes
    • 

  • Genes
    • 

  • Quantization
    • 

  • Population
    • 



This second part introduces the concept and implements of genetic operations (cross-over, mutation and selection) on a population of chromosomes. These operators are applied recursively on each chromosome and each of the genes it contains.

The 3rd and final post on genetic algorithms, explores the application of genetic algorithm as a solver or optimizer Genetic Algorithms III: Solver








Note: 
For the sake of readability of the implementation of algorithms, all non-essential code such as error checking, comments, exception, validation of class and method arguments, scoping qualifiers or import is omitted




Selection






The first genetic operator of the reproduction cycle is the selection process. The
select method of the class Population implements the steps of the selection of the fittest chromosomes in the population in the most efficient manner, as follows:







 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
def select(score: Chromosome[T] => Unit, cutOff: Double) = {
   val cumul = chromosomes./:(0.0)(
      (s,x) =>{ score(xy); s + xy. unfitness} 
   ) 

   chromosomes foreach( _ /= cumul) 
   val newChromosomes = chromosomes.sortWith(_.unfitness < _.unfitness)
   val cutOffSize = (cutOff*newChromosomes.size).floor.toInt
   val newPopSize = if(limit<cutOffSize) limit else cutOffSize

   chromosomes.clear
   chromosomes ++= newChromosomes.take(newPopSize)
}











The select method computes the cumulative sum of an unfitness value, cumul,
for the entire population (lines 2 -3). It normalizes the unfitness of each chromosome (line 6), orders the population by decreasing value (line 7), and applies a soft limit function on population growth, cutOff (line 8). The last step reduces the size of the population to the lowest of the two limits: the hard limit, limit, or the soft limit,cutOffSize (line 9).




Cross-over


There
are several options to select pairs of chromosomes for crossover. This implementation
ranks the chromosomes by their fitness value and then divides the population into two
halves. Finally, it pairs the chromosomes of identical rank from each half as illustrated
in the following diagram: 















Population cross-over

The crossover implementation, +- , selects the parent chromosome candidates for
crossover using the pairing scheme described earlier.







 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
def +- (xOver: Double): Unit = {
 if( size > 1) {
   val mid = size>>1
   val bottom = chromosomes.slice(mid, size) 
   
   val gIdx = geneticIndices(xOver)
   val offSprings = chromosomes.take(mid)
      .zip(bottom)
      .map(p => p._1 +-(p._2, gIdx))
      .unzip

   chromosomes ++= offSprings._1 ++ offSprings._2
 }
}











The implementation of the cross-over on the population of chromosomes ranked by their unfitness consists of 
    


  • Get the mid point of the list of ranked chromosomes (line 3)
    • 

  • Get the least fit half of the chromosome (line 4)
    • 

  • Retrieve the position of the bit in the chromosome the cross-over applies (line 6)
    • 

  • Retrieve the two offspring by crossing over pairs of chromosomes from each half of the ranked population (lines 7 - 10)
    • 

  • Add the two off-springs to the current population (line 12)
    • 











Chromosome cross-over

The implementation of the crossover for a pair of chromosomes using hierarchical
encoding follows two steps:
    

  • Find the gene on each chromosome that corresponds to the crossover index,
gIdx.chOpIdx, and then swap the remaining genes (line 6)
    • 

  • Split and splice the gene crossover at xoverIdx (lines 8 & 12)
    • 









 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
def +-(
   that: Chromosome[T], 
   gIdx: GeneticIndices
): (Chromosome[T], Chromosome[T]) = {
 
    val xoverIdx = gIdx.chOpIdx
    val xGenes = spliceGene(gIdx, that.code(xoverIdx) )
 
    val offSprng1 = code.slice(0, xoverIdx) ::: xGenes._1 
        :: that.code.drop(xoverIdx+1)

    val offSprng2 = that.code.slice(0, xoverIdx) ::: xGenes._2 
        :: code.drop(xoverIdx+1)
 
   (Chromosome[T](offSprng1), Chromosome[T](offSprng2)
}











The crossover method computes the index of the bit that defies the crossover
xoverIdx in each parent chromosome. The genes code(xoverIdx)
and that.code(xoverIdx)  are swapped and spliced by the spliceGene method
to generate a spliced gene (lines 9 - 13).



The method spliceGene is implemented below.







 
def spliceGene(gIdx: GeneticIndices, thatCode: T): (T, T) = {
     ((this.code(gIdx.chOpIdx) +- (thatCode, gIdx)),
      (thatCode +- (code(gIdx.chOpIdx), gIdx)) )
}










Gene cross-over
The crossover is applied to a gene through the +- method of the Gene class. The
exchange of bits between the two genes this and that uses the BitSet Java class
to rearrange the bits after the permutation (lines 4 - 6). The bit string is then decoded to produce the predicate or logical representation of the gene (line 8).







 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
def +- (that: Gene, idx: GeneticIndices): Gene = {
    val clonedBits = cloneBits(bits)
  
    Range(gIdx.geneOpIdx, bits.size).foreach(n =>
        if( that.bits.get(n) ) clonedBits.set(n)
        else clonedBits.clear(n)
    ) 
    val valOp = decode(clonedBits)

    Gene(id, valOp._1, valOp._2)
}








Mutation 





Population mutation
The mutation operator ^ invokes the same operator for all the chromosomes in the
population and then adds the mutated chromosomes to the existing population,
so that they can compete with the original chromosomes. 
The mutation parameter mu is used to compute the absolute index of the mutating
gene, geneticIndices(mu).
We use the notation ^ to
define the mutation operator to remind the reader that the mutation is implemented
by flipping one bit:






 
def ^ (mu: Double): Unit =
     chromosomes ++= chromosomes.map(_ ^ geneticIndices(mu)) 
 









Chromosome mutation
The implementation of the mutation operator ^ on a chromosome consists of
mutating the gene of the index gIdx.chOpIdx  and then updating the list xs of
genes in the chromosome. The method returns a new chromosome with this new generic code
that is added to the population.








def ^ (gIdx: GeneticIndices): Chromosome[T] = {
    val mutated = code(gIdx.chOpIdx) ^ gIdx

    val xs = Range(0, code.size).map(
        i => if(i==gIdx.chOpIdx) mutated else code(i)
    ).toList

    Chromosome[T](xs)
}










Gene mutation
Finally, the mutation operator flips (XOR) the bit at the index gIdx.geneOpIdx


The ^ method mutates the cloned bit string, clonedBits (line 2) by flipping the bit at
the index gIdx.geneOpIdx (line 3). It decodes line 5) and converts the mutated bit string by
converting it into a (target value, operator) tuple (line 7). The last step creates a new
gene from the target-operator tuple.





1
2
3
4
5
6
7
8
def ^ (gIdx: GeneticIndices): Gene = {
     val clonedBits = cloneBits(bits)
     clonedBits.flip(gIdx.geneOpIdx)

     val valOp = decode(clonedBits)

     Gene(id, valOp._1, valOp._2)
}






This concludes the second part of the implementation of genetic algorithms in Scala, dedicated to genetic operators.

Thank you for reading this article. For more information ...




References


---------------------------
Patrick Nicolas has over 25 years of experience in software and data engineering, architecture design and end-to-end deployment and support with extensive knowledge in machine learning. 
He has been director of data engineering at Aideo Technologies since 2017 and he is the author of "Scala for Machine Learning" Packt Publishing ISBN 978-1-78712-238-3










Posted by






No comments:
  






















Labels:
,
,
,
,
,
,
,
,
















        

      









Subscribe to:
Posts (Atom)