gonum.org/v1/gonum@v0.14.0/stat/mds/mds.go (about)

     1  // Copyright ©2018 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package mds
     6  
     7  import (
     8  	"math"
     9  
    10  	"gonum.org/v1/gonum/blas/blas64"
    11  	"gonum.org/v1/gonum/mat"
    12  )
    13  
    14  // TorgersonScaling converts a dissimilarity matrix to a matrix containing
    15  // Euclidean coordinates. TorgersonScaling places the coordinates in dst and
    16  // returns it and the number of positive Eigenvalues if successful.
    17  // Note that Eigen Decomposition is numerically unstable and so Eigenvalues
    18  // near zero should be examined and the value returned for k is advisory only.
    19  // If the scaling is not successful, dst will be empty upon return.
    20  // When the scaling is successful, dst will be resized to k columns wide.
    21  // Eigenvalues will be copied into eigdst and returned as eig if it is provided.
    22  //
    23  // TorgersonScaling will panic if dst is not empty.
    24  func TorgersonScaling(dst *mat.Dense, eigdst []float64, dis mat.Symmetric) (k int, eig []float64) {
    25  	// https://doi.org/10.1007/0-387-28981-X_12
    26  
    27  	n := dis.SymmetricDim()
    28  	if dst.IsEmpty() {
    29  		dst.ReuseAs(n, n)
    30  	} else {
    31  		panic("mds: receiver matrix not empty")
    32  	}
    33  
    34  	b := mat.NewSymDense(n, nil)
    35  	for i := 0; i < n; i++ {
    36  		for j := i; j < n; j++ {
    37  			v := dis.At(i, j)
    38  			v *= v
    39  			b.SetSym(i, j, v)
    40  		}
    41  	}
    42  	c := mat.NewSymDense(n, nil)
    43  	s := -1 / float64(n)
    44  	for i := 0; i < n; i++ {
    45  		c.SetSym(i, i, 1+s)
    46  		for j := i + 1; j < n; j++ {
    47  			c.SetSym(i, j, s)
    48  		}
    49  	}
    50  	dst.Product(c, b, c)
    51  	for i := 0; i < n; i++ {
    52  		for j := i; j < n; j++ {
    53  			b.SetSym(i, j, -0.5*dst.At(i, j))
    54  		}
    55  	}
    56  
    57  	var ed mat.EigenSym
    58  	ok := ed.Factorize(b, true)
    59  	if !ok {
    60  		return 0, eigdst
    61  	}
    62  	ed.VectorsTo(dst)
    63  	vals := ed.Values(nil)
    64  	reverse(vals, dst.RawMatrix())
    65  	copy(eigdst, vals)
    66  
    67  	for i, v := range vals {
    68  		if v < 0 {
    69  			vals[i] = 0
    70  			continue
    71  		}
    72  		k = i + 1
    73  		vals[i] = math.Sqrt(v)
    74  	}
    75  
    76  	var tmp mat.Dense
    77  	tmp.Mul(dst, mat.NewDiagonalRect(n, k, vals[:k]))
    78  	*dst = *dst.Slice(0, n, 0, k).(*mat.Dense)
    79  	dst.Copy(&tmp)
    80  
    81  	return k, eigdst
    82  }
    83  
    84  func reverse(values []float64, vectors blas64.General) {
    85  	for i, j := 0, len(values)-1; i < j; i, j = i+1, j-1 {
    86  		values[i], values[j] = values[j], values[i]
    87  		blas64.Swap(blas64.Vector{N: vectors.Rows, Inc: vectors.Stride, Data: vectors.Data[i:]},
    88  			blas64.Vector{N: vectors.Rows, Inc: vectors.Stride, Data: vectors.Data[j:]})
    89  	}
    90  }