mat64 NMF ( translated from go.matrix NMF https://gist.github.com/r9y9/9030922 )

mat64_nmf.go

package main

import (
        "fmt"
        . "github.com/gonum/matrix/mat64"
)

func NewElemDense(r int, c int, v float64) *Dense {
        data := make([]float64, r*c)
        for i := range data {
                data[i] = v
        }
        a := NewDense(r, c, data)
        return a
}

func NMF(Y *Dense, numBasis int, numIter int) (*Dense, *Dense) {
        row, col := Y.Dims()
        H := NewElemDense(row, numBasis, 1.0)
        U := NewElemDense(numBasis, col, 1.0)
        epsiron1 := NewElemDense(numBasis, col, 1.0e-07)
        epsiron2 := NewElemDense(row, numBasis, 1.0e-07)
        t2 := NewDense(row, numBasis, nil)
        numerH := NewDense(row, numBasis, nil)
        denomH := NewDense(row, numBasis, nil)
        denomH0 := NewDense(row, numBasis, nil)
        u2 := NewDense(numBasis, numBasis, nil)
        t1 := NewDense(numBasis, col, nil)
        numerU := NewDense(numBasis, col, nil)
        denomU := NewDense(numBasis, col, nil)
        denomU0 := NewDense(numBasis, col, nil)
        h2 := NewDense(numBasis, numBasis, nil)
        product := NewDense(row, col, nil)
        r := NewDense(row, col, nil)
        for iter := 0; iter < numIter; iter++ {
                
                product.Mul(H, U)
                r.Sub(Y, product)
                rv := Norm(r, 2)
                fmt.Println("objective function: ", rv)
                // Update U
                t1.Mul(H.T(), Y)
                numerU.MulElem(U, t1)
                h2.Mul(H.T(), H)
                denomU0.Mul(h2, U)
                denomU.Add(denomU0, epsiron1)
                U.DivElem(numerU, denomU)
                                // Update H
                t2.Mul(Y, U.T())
                numerH.MulElem(H, t2)
                u2.Mul(U, U.T())
                denomH0.Mul(H, u2)
                denomH.Add(denomH0, epsiron2)
                H.DivElem(numerH, denomH)

        }
        return H, U
}

func main() {
        A := NewElemDense(5, 5, 50.0)
        fmt.Println("The original matrix: \n", A)
        H, U := NMF(A, 2, 10)
        fmt.Println("Basis: \n", H)
        fmt.Println("Activity: \n", U)
}