blob: 64dab97d318348915bc477303854e00b0f93c3aa [file] [log] [blame]
 // Copyright ©2013 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package mat import ( "testing" "golang.org/x/exp/rand" ) func TestLQ(t *testing.T) { t.Parallel() rnd := rand.New(rand.NewSource(1)) for _, test := range []struct { m, n int }{ {5, 5}, {5, 10}, } { m := test.m n := test.n a := NewDense(m, n, nil) for i := 0; i < m; i++ { for j := 0; j < n; j++ { a.Set(i, j, rnd.NormFloat64()) } } var want Dense want.CloneFrom(a) var lq LQ lq.Factorize(a) var l, q Dense lq.QTo(&q) if !isOrthonormal(&q, 1e-10) { t.Errorf("Q is not orthonormal: m = %v, n = %v", m, n) } lq.LTo(&l) var got Dense got.Mul(&l, &q) if !EqualApprox(&got, &want, 1e-12) { t.Errorf("LQ does not equal original matrix. \nWant: %v\nGot: %v", want, got) } } } func TestLQSolveTo(t *testing.T) { t.Parallel() rnd := rand.New(rand.NewSource(1)) for _, trans := range []bool{false, true} { for _, test := range []struct { m, n, bc int }{ {5, 5, 1}, {5, 10, 1}, {5, 5, 3}, {5, 10, 3}, } { m := test.m n := test.n bc := test.bc a := NewDense(m, n, nil) for i := 0; i < m; i++ { for j := 0; j < n; j++ { a.Set(i, j, rnd.Float64()) } } br := m if trans { br = n } b := NewDense(br, bc, nil) for i := 0; i < br; i++ { for j := 0; j < bc; j++ { b.Set(i, j, rnd.Float64()) } } var x Dense lq := &LQ{} lq.Factorize(a) err := lq.SolveTo(&x, trans, b) if err != nil { t.Errorf("unexpected error from LQ solve: %v", err) } // Test that the normal equations hold. // Aᵀ * A * x = Aᵀ * b if !trans // A * Aᵀ * x = A * b if trans var lhs Dense var rhs Dense if trans { var tmp Dense tmp.Mul(a, a.T()) lhs.Mul(&tmp, &x) rhs.Mul(a, b) } else { var tmp Dense tmp.Mul(a.T(), a) lhs.Mul(&tmp, &x) rhs.Mul(a.T(), b) } if !EqualApprox(&lhs, &rhs, 1e-10) { t.Errorf("Normal equations do not hold.\nLHS: %v\n, RHS: %v\n", lhs, rhs) } } } // TODO(btracey): Add in testOneInput when it exists. } func TestLQSolveToVec(t *testing.T) { t.Parallel() rnd := rand.New(rand.NewSource(1)) for _, trans := range []bool{false, true} { for _, test := range []struct { m, n int }{ {5, 5}, {5, 10}, } { m := test.m n := test.n a := NewDense(m, n, nil) for i := 0; i < m; i++ { for j := 0; j < n; j++ { a.Set(i, j, rnd.Float64()) } } br := m if trans { br = n } b := NewVecDense(br, nil) for i := 0; i < br; i++ { b.SetVec(i, rnd.Float64()) } var x VecDense lq := &LQ{} lq.Factorize(a) err := lq.SolveVecTo(&x, trans, b) if err != nil { t.Errorf("unexpected error from LQ solve: %v", err) } // Test that the normal equations hold. // Aᵀ * A * x = Aᵀ * b if !trans // A * Aᵀ * x = A * b if trans var lhs Dense var rhs Dense if trans { var tmp Dense tmp.Mul(a, a.T()) lhs.Mul(&tmp, &x) rhs.Mul(a, b) } else { var tmp Dense tmp.Mul(a.T(), a) lhs.Mul(&tmp, &x) rhs.Mul(a.T(), b) } if !EqualApprox(&lhs, &rhs, 1e-10) { t.Errorf("Normal equations do not hold.\nLHS: %v\n, RHS: %v\n", lhs, rhs) } } } // TODO(btracey): Add in testOneInput when it exists. } func TestLQSolveToCond(t *testing.T) { t.Parallel() for _, test := range []*Dense{ NewDense(2, 2, []float64{1, 0, 0, 1e-20}), NewDense(2, 3, []float64{1, 0, 0, 0, 1e-20, 0}), } { m, _ := test.Dims() var lq LQ lq.Factorize(test) b := NewDense(m, 2, nil) var x Dense if err := lq.SolveTo(&x, false, b); err == nil { t.Error("No error for near-singular matrix in matrix solve.") } bvec := NewVecDense(m, nil) var xvec VecDense if err := lq.SolveVecTo(&xvec, false, bvec); err == nil { t.Error("No error for near-singular matrix in matrix solve.") } } }