calculating least squares

foo.go

package main

import (
	"code.google.com/p/gomatrix/matrix"
	"fmt"
	"log"
)

func solve(x, b *matrix.DenseMatrix) (y *matrix.DenseMatrix, err error) {
	if x.Cols() > x.Rows() {
		x = x.Transpose()
		defer func() {
			y = y.Transpose()
		}()
	}
	u, sig, v, err := x.SVD()
	if err != nil {
		return
	}

	//invert the diagonal
	for i := 0; i < sig.Cols(); i++ {
		sig.Set(i, i, 1.0/sig.Get(i, i))
	}
	//transpose
	if err = sig.TransposeInPlace(); err != nil {
		return
	}

	s1, err := u.Transpose().TimesDense(b)
	if err != nil {
		return
	}
	
	//  Diagonal multiplication in-place
	for i := 0; i < sig.Cols(); i++ {
		s1.Set(i, 0, sig.Get(i, i) * s1.Get(i, 0))
        }

	y, err = v.TimesDense(s1)
	if err != nil {
		return
	}
	return
}

func main() {
	a := matrix.MakeDenseMatrixStacked([][]float64{
		{1, 0, 1},
		{0, 1, 1},
		{0, 0, 1},
		{1, 0, 1},
		{0, 1, 1},
		{0, 0, 1},
	})
	b := matrix.MakeDenseMatrix([]float64{
		14,
		20,
		3,
		14,
		19,
		4,
	}, 6, 1)

	y, err := solve(a, b)
	if err != nil {
		log.Fatal(err)
	}
	fmt.Println(y)
}