mat: Change Eigen to use complex matrix representation
* mat: Change Eigen to use complex matrix representation
Fixes
MODIFY
mat/cmatrix.go
mat/cmatrix.go
package mat
-import "math/cmplx"
+import (
+ "math"
+ "math/cmplx"
+
+ "gonum.org/v1/gonum/floats"
+)
// CMatrix is the basic matrix interface type for complex matrices.
type CMatrix interface {
@@ -71,3 +76,135 @@ type Unconjugator interface {
// conjugate transpose.
Unconjugate() CMatrix
}
+
+// useC returns a complex128 slice with l elements, using c if it
+// has the necessary capacity, otherwise creating a new slice.
+func useC(c []complex128, l int) []complex128 {
+ if l <= cap(c) {
+ return c[:l]
+ }
+ return make([]complex128, l)
+}
+
+// useZeroedC returns a complex128 slice with l elements, using c if it
+// has the necessary capacity, otherwise creating a new slice. The
+// elements of the returned slice are guaranteed to be zero.
+func useZeroedC(c []complex128, l int) []complex128 {
+ if l <= cap(c) {
+ c = c[:l]
+ zeroC(c)
+ return c
+ }
+ return make([]complex128, l)
+}
+
+// zeroC zeros the given slice's elements.
+func zeroC(c []complex128) {
+ for i := range c {
+ c[i] = 0
+ }
+}
+
+// unconjugate unconjugates a matrix if applicable. If a is an Unconjugator, then
+// unconjugate returns the underlying matrix and true. If it is not, then it returns
+// the input matrix and false.
+func unconjugate(a CMatrix) (CMatrix, bool) {
+ if ut, ok := a.(Unconjugator); ok {
+ return ut.Unconjugate(), true
+ }
+ return a, false
+}
+
+// CEqual returns whether the matrices a and b have the same size
+// and are element-wise equal.
+func CEqual(a, b CMatrix) bool {
+ ar, ac := a.Dims()
+ br, bc := b.Dims()
+ if ar != br || ac != bc {
+ return false
+ }
+ // TODO(btracey): Add in fast-paths.
+ for i := 0; i < ar; i++ {
+ for j := 0; j < ac; j++ {
+ if a.At(i, j) != b.At(i, j) {
+ return false
+ }
+ }
+ }
+ return true
+}
+
+// CEqualApprox returns whether the matrices a and b have the same size and contain all equal
+// elements with tolerance for element-wise equality specified by epsilon. Matrices
+// with non-equal shapes are not equal.
+func CEqualApprox(a, b CMatrix, epsilon float64) bool {
+ // TODO(btracey):
+ ar, ac := a.Dims()
+ br, bc := b.Dims()
+ if ar != br || ac != bc {
+ return false
+ }
+ for i := 0; i < ar; i++ {
+ for j := 0; j < ac; j++ {
+ if !cEqualWithinAbsOrRel(a.At(i, j), b.At(i, j), epsilon, epsilon) {
+ return false
+ }
+ }
+ }
+ return true
+}
+
+// TODO(btracey): Move these into a cmplxs if/when we have one.
+
+func cEqualWithinAbsOrRel(a, b complex128, absTol, relTol float64) bool {
+ if cEqualWithinAbs(a, b, absTol) {
+ return true
+ }
+ return cEqualWithinRel(a, b, relTol)
+}
+
+// cEqualWithinAbs returns true if a and b have an absolute
+// difference of less than tol.
+func cEqualWithinAbs(a, b complex128, tol float64) bool {
+ return a == b || cmplx.Abs(a-b) <= tol
+}
+
+const minNormalFloat64 = 2.2250738585072014e-308
+
+// cEqualWithinRel returns true if the difference between a and b
+// is not greater than tol times the greater value.
+func cEqualWithinRel(a, b complex128, tol float64) bool {
+ if a == b {
+ return true
+ }
+ if cmplx.IsNaN(a) || cmplx.IsNaN(b) {
+ return false
+ }
+ // Cannot play the same trick as in floats because there are multiple
+ // possible infinities.
+ if cmplx.IsInf(a) {
+ if !cmplx.IsInf(b) {
+ return false
+ }
+ ra := real(a)
+ if math.IsInf(ra, 0) {
+ if ra == real(b) {
+ return floats.EqualWithinRel(imag(a), imag(b), tol)
+ }
+ return false
+ }
+ if imag(a) == imag(b) {
+ return floats.EqualWithinRel(ra, real(b), tol)
+ }
+ return false
+ }
+ if cmplx.IsInf(b) {
+ return false
+ }
+
+ delta := cmplx.Abs(a - b)
+ if delta <= minNormalFloat64 {
+ return delta <= tol*minNormalFloat64
+ }
+ return delta/math.Max(cmplx.Abs(a), cmplx.Abs(b)) <= tol
+}

MODIFY
mat/eigen_test.go
mat/eigen_test.go
@@ -18,8 +18,8 @@ func TestEigen(t *testing.T) {
a *Dense
values []complex128
- left *Dense
- right *Dense
+ left *CDense
+ right *CDense
}{
{
a: NewDense(3, 3, []float64{
@@ -28,17 +28,39 @@ func TestEigen(t *testing.T) {
0, 0, 1,
}),
values: []complex128{1, 1, 1},
- left: NewDense(3, 3, []float64{
+ left: NewCDense(3, 3, []complex128{
1, 0, 0,
0, 1, 0,
0, 0, 1,
}),
- right: NewDense(3, 3, []float64{
+ right: NewCDense(3, 3, []complex128{
1, 0, 0,
0, 1, 0,
0, 0, 1,
}),
},
+ {
+ // Values compared with numpy.
+ a: NewDense(4, 4, []float64{
+ 0.9025, 0.025, 0.475, 0.0475,
+ 0.0475, 0.475, 0.475, 0.0025,
+ 0.0475, 0.025, 0.025, 0.9025,
+ 0.0025, 0.475, 0.025, 0.0475,
+ }),
+ values: []complex128{1, 0.7300317046114154, -0.1400158523057075 + 0.452854925738716i, -0.1400158523057075 - 0.452854925738716i},
+ left: NewCDense(4, 4, []complex128{
+ 0.5, -0.3135167160788314, 0.0205812178013689 - 0.0045809393001271i, 0.0205812178013689 + 0.0045809393001271i,
+ 0.5, 0.7842199280224774, -0.3755102695419336 + 0.2924634904103882i, -0.3755102695419336 - 0.2924634904103882i,
+ 0.5, 0.3320220078078358, -0.1605261632278496 - 0.3881393645202528i, -0.1605261632278496 + 0.3881393645202528i,
+ 0.5, 0.4200806584012395, 0.7723935249234153, 0.7723935249234153,
+ }),
+ right: NewCDense(4, 4, []complex128{
+ 0.9476399565969628, -0.8637347682162745, -0.2688989440320280 - 0.1282234938321029i, -0.2688989440320280 + 0.1282234938321029i,
+ 0.2394935907064427, 0.3457075153704627, -0.3621360383713332 - 0.2583198964498771i, -0.3621360383713332 + 0.2583198964498771i,
+ 0.1692743801716332, 0.2706851011641580, 0.7426369401030960, 0.7426369401030960,
+ 0.1263626404003607, 0.2473421516816520, -0.1116019576997347 + 0.3865433902819795i, -0.1116019576997347 - 0.3865433902819795i,
+ }),
+ },
} {
var e1, e2, e3, e4 Eigen
ok := e1.Factorize(test.a, true, true)
@@ -50,13 +72,13 @@ func TestEigen(t *testing.T) {
e4.Factorize(test.a, false, false)
v1 := e1.Values(nil)
- if !cmplxEqual(v1, test.values) {
- t.Errorf("eigenvector mismatch. Case %v", i)
+ if !cmplxEqualTol(v1, test.values, 1e-14) {
+ t.Errorf("eigenvalue mismatch. Case %v", i)
}
- if !Equal(e1.LeftVectors(), test.left) {
+ if !CEqualApprox(e1.LeftVectorsTo(nil), test.left, 1e-14) {
t.Errorf("left eigenvector mismatch. Case %v", i)
}
- if !Equal(e1.Vectors(), test.right) {
+ if !CEqualApprox(e1.VectorsTo(nil), test.right, 1e-14) {
t.Errorf("right eigenvector mismatch. Case %v", i)
}
@@ -70,11 +92,11 @@ func TestEigen(t *testing.T) {
if !cmplxEqual(v1, e4.Values(nil)) {
t.Errorf("eigenvector mismatch. Case %v", i)
}
- if !Equal(e1.Vectors(), e2.Vectors()) {
+ if !CEqual(e1.VectorsTo(nil), e2.VectorsTo(nil)) {
t.Errorf("right eigenvector mismatch. Case %v", i)
}
- if !Equal(e1.LeftVectors(), e3.LeftVectors()) {
- t.Errorf("right eigenvector mismatch. Case %v", i)
+ if !CEqual(e1.LeftVectorsTo(nil), e3.LeftVectorsTo(nil)) {
+ t.Errorf("left eigenvector mismatch. Case %v", i)
}
// TODO(btracey): Also add in a test for correctness when #308 is
@@ -91,6 +113,15 @@ func cmplxEqual(v1, v2 []complex128) bool {
return true
}
+func cmplxEqualTol(v1, v2 []complex128, tol float64) bool {
+ for i, v := range v1 {
+ if !cEqualWithinAbsOrRel(v, v2[i], tol, tol) {
+ return false
+ }
+ }
+ return true
+}
+
func TestSymEigen(t *testing.T) {
// Hand coded tests with results from lapack.
for _, test := range []struct {

MODIFY
mat/hogsvd.go
mat/hogsvd.go
@@ -101,7 +101,18 @@ func (gsvd *HOGSVD) Factorize(m ...Matrix) (ok bool) {
gsvd.err = errors.New("hogsvd: eigen decomposition failed")
return false
}
- v := eig.Vectors()
+ vc := eig.VectorsTo(nil)
+ // vc is guaranteed to have real eigenvalues.
+ rc, cc := vc.Dims()
+ v := NewDense(rc, cc, nil)
+ for i := 0; i < rc; i++ {
+ for j := 0; j < cc; j++ {
+ a := vc.At(i, j)
+ v.set(i, j, real(a))
+ }
+ }
+ // Rescale the columns of v by their Frobenius norms.
+ // Work done in cv is reflected in v.
var cv VecDense
for j := 0; j < c; j++ {
cv.ColViewOf(v, j)

mat: Change Eigen to use complex matrix representation (#849)

* mat: Change Eigen to use complex matrix representation

Fixes #738