blob: cc047b478b95b511b89bc18681910d2019824daf [file] [log] [blame]
 // Copyright ©2015 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 testlapack import ( "fmt" "math" "sort" "testing" "golang.org/x/exp/rand" "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/blas/blas64" "gonum.org/v1/gonum/floats" "gonum.org/v1/gonum/lapack" ) type Dgesvder interface { Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float64, lda int, s, u []float64, ldu int, vt []float64, ldvt int, work []float64, lwork int) (ok bool) } func DgesvdTest(t *testing.T, impl Dgesvder, tol float64) { for _, m := range []int{0, 1, 2, 3, 4, 5, 10, 150, 300} { for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 150} { for _, mtype := range []int{1, 2, 3, 4, 5} { dgesvdTest(t, impl, m, n, mtype, tol) } } } } // dgesvdTest tests a Dgesvd implementation on an m×n matrix A generated // according to mtype as: // - the zero matrix if mtype == 1, // - the identity matrix if mtype == 2, // - a random matrix with a given condition number and singular values if mtype == 3, 4, or 5. // It first computes the full SVD A = U*Sigma*Vᵀ and checks that // - U has orthonormal columns, and Vᵀ has orthonormal rows, // - U*Sigma*Vᵀ multiply back to A, // - the singular values are non-negative and sorted in decreasing order. // Then all combinations of partial SVD results are computed and checked whether // they match the full SVD result. func dgesvdTest(t *testing.T, impl Dgesvder, m, n, mtype int, tol float64) { const tolOrtho = 1e-15 rnd := rand.New(rand.NewSource(1)) // Use a fixed leading dimension to reduce testing time. lda := n + 3 ldu := m + 5 ldvt := n + 7 minmn := min(m, n) // Allocate A and fill it with random values. The in-range elements will // be overwritten below according to mtype. a := make([]float64, m*lda) for i := range a { a[i] = rnd.NormFloat64() } var aNorm float64 switch mtype { default: panic("unknown test matrix type") case 1: // Zero matrix. for i := 0; i < m; i++ { for j := 0; j < n; j++ { a[i*lda+j] = 0 } } aNorm = 0 case 2: // Identity matrix. for i := 0; i < m; i++ { for j := 0; j < n; j++ { if i == j { a[i*lda+i] = 1 } else { a[i*lda+j] = 0 } } } aNorm = 1 case 3, 4, 5: // Scaled random matrix. // Generate singular values. s := make([]float64, minmn) Dlatm1(s, 4, // s[i] = 1 - i*(1-1/cond)/(minmn-1) float64(max(1, minmn)), // where cond = max(1,minmn) false, // signs of s[i] are not randomly flipped 1, rnd) // random numbers are drawn uniformly from [0,1) // Decide scale factor for the singular values based on the matrix type. ulp := dlamchP unfl := dlamchS ovfl := 1 / unfl aNorm = 1 if mtype == 4 { aNorm = unfl / ulp } if mtype == 5 { aNorm = ovfl * ulp } // Scale singular values so that the maximum singular value is // equal to aNorm (we know that the singular values are // generated above to be spread linearly between 1/cond and 1). floats.Scale(aNorm, s) // Generate A by multiplying S by random orthogonal matrices // from left and right. Dlagge(m, n, max(0, m-1), max(0, n-1), s, a, lda, rnd, make([]float64, m+n)) } aCopy := make([]float64, len(a)) copy(aCopy, a) for _, wl := range []worklen{minimumWork, mediumWork, optimumWork} { // Restore A because Dgesvd overwrites it. copy(a, aCopy) // Allocate slices that will be used below to store the results of full // SVD and fill them. uAll := make([]float64, m*ldu) for i := range uAll { uAll[i] = rnd.NormFloat64() } vtAll := make([]float64, n*ldvt) for i := range vtAll { vtAll[i] = rnd.NormFloat64() } sAll := make([]float64, min(m, n)) for i := range sAll { sAll[i] = math.NaN() } prefix := fmt.Sprintf("m=%v,n=%v,work=%v,mtype=%v", m, n, wl, mtype) // Determine workspace size based on wl. minwork := max(1, max(5*min(m, n), 3*min(m, n)+max(m, n))) var lwork int switch wl { case minimumWork: lwork = minwork case mediumWork: work := make([]float64, 1) impl.Dgesvd(lapack.SVDAll, lapack.SVDAll, m, n, a, lda, sAll, uAll, ldu, vtAll, ldvt, work, -1) lwork = (int(work[0]) + minwork) / 2 case optimumWork: work := make([]float64, 1) impl.Dgesvd(lapack.SVDAll, lapack.SVDAll, m, n, a, lda, sAll, uAll, ldu, vtAll, ldvt, work, -1) lwork = int(work[0]) } work := make([]float64, max(1, lwork)) for i := range work { work[i] = math.NaN() } // Compute the full SVD which will be used later for checking the partial results. ok := impl.Dgesvd(lapack.SVDAll, lapack.SVDAll, m, n, a, lda, sAll, uAll, ldu, vtAll, ldvt, work, len(work)) if !ok { t.Fatalf("Case %v: unexpected failure in full SVD", prefix) } // Check that uAll, sAll, and vtAll multiply back to A by computing a residual // |A - U*S*VT| / (n*aNorm) if resid := svdFullResidual(m, n, aNorm, aCopy, lda, uAll, ldu, sAll, vtAll, ldvt); resid > tol { t.Errorf("Case %v: original matrix not recovered for full SVD, |A - U*D*VT|=%v", prefix, resid) } if minmn > 0 { // Check that uAll is orthogonal. q := blas64.General{Rows: m, Cols: m, Data: uAll, Stride: ldu} if resid := residualOrthogonal(q, false); resid > tolOrtho*float64(m) { t.Errorf("Case %v: UAll is not orthogonal; resid=%v, want<=%v", prefix, resid, tolOrtho*float64(m)) } // Check that vtAll is orthogonal. q = blas64.General{Rows: n, Cols: n, Data: vtAll, Stride: ldvt} if resid := residualOrthogonal(q, false); resid > tolOrtho*float64(n) { t.Errorf("Case %v: VTAll is not orthogonal; resid=%v, want<=%v", prefix, resid, tolOrtho*float64(n)) } } // Check that singular values are decreasing. if !sort.IsSorted(sort.Reverse(sort.Float64Slice(sAll))) { t.Errorf("Case %v: singular values from full SVD are not decreasing", prefix) } // Check that singular values are non-negative. if minmn > 0 && floats.Min(sAll) < 0 { t.Errorf("Case %v: some singular values from full SVD are negative", prefix) } // Do partial SVD and compare the results to sAll, uAll, and vtAll. for _, jobU := range []lapack.SVDJob{lapack.SVDAll, lapack.SVDStore, lapack.SVDOverwrite, lapack.SVDNone} { for _, jobVT := range []lapack.SVDJob{lapack.SVDAll, lapack.SVDStore, lapack.SVDOverwrite, lapack.SVDNone} { if jobU == lapack.SVDOverwrite || jobVT == lapack.SVDOverwrite { // Not implemented. continue } if jobU == lapack.SVDAll && jobVT == lapack.SVDAll { // Already checked above. continue } prefix := prefix + ",job=" + svdJobString(jobU) + "U-" + svdJobString(jobVT) + "VT" // Restore A to its original values. copy(a, aCopy) // Allocate slices for the results of partial SVD and fill them. u := make([]float64, m*ldu) for i := range u { u[i] = rnd.NormFloat64() } vt := make([]float64, n*ldvt) for i := range vt { vt[i] = rnd.NormFloat64() } s := make([]float64, min(m, n)) for i := range s { s[i] = math.NaN() } for i := range work { work[i] = math.NaN() } ok := impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, u, ldu, vt, ldvt, work, len(work)) if !ok { t.Fatalf("Case %v: unexpected failure in partial Dgesvd", prefix) } if minmn == 0 { // No panic and the result is ok, there is // nothing else to check. continue } // Check that U has orthogonal columns and that it matches UAll. switch jobU { case lapack.SVDStore: q := blas64.General{Rows: m, Cols: minmn, Data: u, Stride: ldu} if resid := residualOrthogonal(q, false); resid > tolOrtho*float64(m) { t.Errorf("Case %v: columns of U are not orthogonal; resid=%v, want<=%v", prefix, resid, tolOrtho*float64(m)) } if res := svdPartialUResidual(m, minmn, u, uAll, ldu); res > tol { t.Errorf("Case %v: columns of U do not match UAll", prefix) } case lapack.SVDAll: q := blas64.General{Rows: m, Cols: m, Data: u, Stride: ldu} if resid := residualOrthogonal(q, false); resid > tolOrtho*float64(m) { t.Errorf("Case %v: columns of U are not orthogonal; resid=%v, want<=%v", prefix, resid, tolOrtho*float64(m)) } if res := svdPartialUResidual(m, m, u, uAll, ldu); res > tol { t.Errorf("Case %v: columns of U do not match UAll", prefix) } } // Check that VT has orthogonal rows and that it matches VTAll. switch jobVT { case lapack.SVDStore: q := blas64.General{Rows: minmn, Cols: n, Data: vtAll, Stride: ldvt} if resid := residualOrthogonal(q, true); resid > tolOrtho*float64(n) { t.Errorf("Case %v: rows of VT are not orthogonal; resid=%v, want<=%v", prefix, resid, tolOrtho*float64(n)) } if res := svdPartialVTResidual(minmn, n, vt, vtAll, ldvt); res > tol { t.Errorf("Case %v: rows of VT do not match VTAll", prefix) } case lapack.SVDAll: q := blas64.General{Rows: n, Cols: n, Data: vtAll, Stride: ldvt} if resid := residualOrthogonal(q, true); resid > tolOrtho*float64(n) { t.Errorf("Case %v: rows of VT are not orthogonal; resid=%v, want<=%v", prefix, resid, tolOrtho*float64(n)) } if res := svdPartialVTResidual(n, n, vt, vtAll, ldvt); res > tol { t.Errorf("Case %v: rows of VT do not match VTAll", prefix) } } // Check that singular values are decreasing. if !sort.IsSorted(sort.Reverse(sort.Float64Slice(s))) { t.Errorf("Case %v: singular values from full SVD are not decreasing", prefix) } // Check that singular values are non-negative. if floats.Min(s) < 0 { t.Errorf("Case %v: some singular values from full SVD are negative", prefix) } if !floats.EqualApprox(s, sAll, tol/10) { t.Errorf("Case %v: singular values differ between full and partial SVD\n%v\n%v", prefix, s, sAll) } } } } } // svdFullResidual returns // |A - U*D*VT| / (n * aNorm) // where U, D, and VT are as computed by Dgesvd with jobU = jobVT = lapack.SVDAll. func svdFullResidual(m, n int, aNorm float64, a []float64, lda int, u []float64, ldu int, d []float64, vt []float64, ldvt int) float64 { // The implementation follows TESTING/dbdt01.f from the reference. minmn := min(m, n) if minmn == 0 { return 0 } // j-th column of A - U*D*VT. aMinusUDVT := make([]float64, m) // D times the j-th column of VT. dvt := make([]float64, minmn) // Compute the residual |A - U*D*VT| one column at a time. var resid float64 for j := 0; j < n; j++ { // Copy j-th column of A to aj. blas64.Copy(blas64.Vector{N: m, Data: a[j:], Inc: lda}, blas64.Vector{N: m, Data: aMinusUDVT, Inc: 1}) // Multiply D times j-th column of VT. for i := 0; i < minmn; i++ { dvt[i] = d[i] * vt[i*ldvt+j] } // Compute the j-th column of A - U*D*VT. blas64.Gemv(blas.NoTrans, -1, blas64.General{Rows: m, Cols: minmn, Data: u, Stride: ldu}, blas64.Vector{N: minmn, Data: dvt, Inc: 1}, 1, blas64.Vector{N: m, Data: aMinusUDVT, Inc: 1}) resid = math.Max(resid, blas64.Asum(blas64.Vector{N: m, Data: aMinusUDVT, Inc: 1})) } if aNorm == 0 { if resid != 0 { // Original matrix A is zero but the residual is non-zero, // return infinity. return math.Inf(1) } // Original matrix A is zero, residual is zero, return 0. return 0 } // Original matrix A is non-zero. if aNorm >= resid { resid = resid / aNorm / float64(n) } else { if aNorm < 1 { resid = math.Min(resid, float64(n)*aNorm) / aNorm / float64(n) } else { resid = math.Min(resid/aNorm, float64(n)) / float64(n) } } return resid } // svdPartialUResidual compares U and URef to see if their columns span the same // spaces. It returns the maximum over columns of // |URef(i) - S*U(i)| // where URef(i) and U(i) are the i-th columns of URef and U, respectively, and // S is ±1 chosen to minimize the expression. func svdPartialUResidual(m, n int, u, uRef []float64, ldu int) float64 { var res float64 for j := 0; j < n; j++ { imax := blas64.Iamax(blas64.Vector{N: m, Data: uRef[j:], Inc: ldu}) s := math.Copysign(1, uRef[imax*ldu+j]) * math.Copysign(1, u[imax*ldu+j]) for i := 0; i < m; i++ { diff := math.Abs(uRef[i*ldu+j] - s*u[i*ldu+j]) res = math.Max(res, diff) } } return res } // svdPartialVTResidual compares VT and VTRef to see if their rows span the same // spaces. It returns the maximum over rows of // |VTRef(i) - S*VT(i)| // where VTRef(i) and VT(i) are the i-th columns of VTRef and VT, respectively, and // S is ±1 chosen to minimize the expression. func svdPartialVTResidual(m, n int, vt, vtRef []float64, ldvt int) float64 { var res float64 for i := 0; i < m; i++ { jmax := blas64.Iamax(blas64.Vector{N: n, Data: vtRef[i*ldvt:], Inc: 1}) s := math.Copysign(1, vtRef[i*ldvt+jmax]) * math.Copysign(1, vt[i*ldvt+jmax]) for j := 0; j < n; j++ { diff := math.Abs(vtRef[i*ldvt+j] - s*vt[i*ldvt+j]) res = math.Max(res, diff) } } return res }