Stability of SVD
Purpose
To explore stability for SVD
> library(Matrix) > N <- 2500 > m <- 50 > n <- 50 > set.seed(12345) > B <- matrix(rnorm(N), nrow = m, ncol = n) > U <- qr.Q(qr(B)) > B <- matrix(rnorm(N), nrow = m, ncol = n) > V <- qr.Q(qr(B)) > d <- diag(sort(runif(50), decreasing = T)) > A <- U %*% d %*% t(V) > U1 <- svd(A)$u > d1 <- diag(svd(A)$d) > V1 <- svd(A)$v |
> normv <- function(x) max(svd(x)$d) |
Norm Comparision for U, V, D
> error.U1 <- normv(U1 - U) > error.V1 <- normv(V1 - V) > error.d1 <- normv(d1 - d) > error <- normv(U1 %*% d1 %*% t(V1) - A) > print(c(error.U1, error.V1, error.d1, error)) [1] 2.000000e+00 2.000000e+00 5.551115e-16 6.463822e-15 |
Fix the errors
> us <- which(abs(U1[1, ] - U[1, ]) > 10^-4) > vs <- which(abs(V1[1, ] - V[1, ]) > 10^-4) > j <- 1 > for (j in seq_along(us)) U1[, us[j]] <- U1[, us[j]] * -1 > for (j in seq_along(vs)) V1[, vs[j]] <- V1[, vs[j]] * -1 |
Norm Comparision for U, V, D
> error.U1.f <- normv(U1 - U) > error.V1.f <- normv(V1 - V) > error.d1.f <- normv(d1 - d) > error.f <- normv(U1 %*% d1 %*% t(V1) - A) > print(c(error.U1.f, error.V1.f, error.d1.f, error.f)) [1] 1.633826e-12 1.633677e-12 5.551115e-16 6.463822e-15 |
- The errors in U1 and V1 are of the order 10^(-12). This is much better than QR algo where QR had a huge error .
Let me write a funciton which does this
> svdtest <- function(input.seed) {
+ set.seed(input.seed)
+ N <- 2500
+ m <- 50
+ n <- 50
+ B <- matrix(rnorm(N), nrow = m, ncol = n)
+ U <- qr.Q(qr(B))
+ B <- matrix(rnorm(N), nrow = m, ncol = n)
+ V <- qr.Q(qr(B))
+ d <- diag(sort(runif(50), decreasing = T))
+ A <- U %*% d %*% t(V)
+ U1 <- svd(A)$u
+ d1 <- diag(svd(A)$d)
+ V1 <- svd(A)$v
+ normv <- function(x) max(svd(x)$d)
+ error.U1 <- normv(U1 - U)
+ error.V1 <- normv(V1 - V)
+ error.d1 <- normv(d1 - d)
+ error <- normv(U1 %*% d1 %*% t(V1) - A)
+ us <- which(abs(U1[1, ] - U[1, ]) > 10^-4)
+ vs <- which(abs(V1[1, ] - V[1, ]) > 10^-4)
+ j <- 1
+ for (j in seq_along(us)) U1[, us[j]] <- U1[, us[j]] * -1
+ for (j in seq_along(vs)) V1[, vs[j]] <- V1[, vs[j]] * -1
+ error.U1.f <- normv(U1 - U)
+ error.V1.f <- normv(V1 - V)
+ error.d1.f <- normv(d1 - d)
+ error.f <- normv(U1 %*% d1 %*% t(V1) - A)
+ res1 <- (c(error.U1, error.V1, error.d1, error))
+ res2 <- (c(error.U1.f, error.V1.f, error.d1.f, error.f))
+ res <- (rbind(res1, res2))
+ cond <- max(diag(d1))/min(diag(d1))
+ result <- list(res = res, cond = cond)
+ return(result)
+ }
> x <- svdtest(123456)
> print(x$res)
[,1] [,2] [,3] [,4]
res1 2.000000e+00 2.000000e+00 5.551115e-16 2.031852e-15
res2 1.108035e-13 1.107732e-13 5.551115e-16 2.031852e-15
> print(x$cond)
[1] 21.08140
> x <- svdtest(12345)
> print(x$res)
[,1] [,2] [,3] [,4]
res1 2.000000e+00 2.000000e+00 5.551115e-16 6.463822e-15
res2 1.633826e-12 1.633677e-12 5.551115e-16 6.463822e-15
> print(x$cond)
[1] 78.51633
> x <- svdtest(1234)
> print(x$res)
[,1] [,2] [,3] [,4]
res1 2.000000e+00 2.000000e+00 1.221245e-15 3.206594e-15
res2 1.891009e-13 1.892900e-13 1.221245e-15 3.206594e-15
> print(x$cond)
[1] 90.96732
> x <- svdtest(123)
> print(x$res)
[,1] [,2] [,3] [,4]
res1 2.000000e+00 2.000000e+00 1.554312e-15 3.897969e-15
res2 8.446307e-13 8.455416e-13 1.554312e-15 3.897969e-15
> print(x$cond)
[1] 44.97811
> x <- svdtest(12)
> print(x$res)
[,1] [,2] [,3] [,4]
res1 2.000000e+00 2.000000e+00 1.221245e-15 3.479001e-15
res2 4.420193e-13 4.421101e-13 1.221245e-15 3.479001e-15
> print(x$cond)
[1] 61.13586 |
Empirically what is observed is for a higher condition number, the errors in U and V are lesser as compared to a lower condition number
Lower condition number means higher eigen values and hence the norm of the error is higher.
Repeat S by raising to 6th power
> svdtest1 <- function(input.seed) {
+ set.seed(input.seed)
+ N <- 2500
+ m <- 50
+ n <- 50
+ B <- matrix(rnorm(N), nrow = m, ncol = n)
+ U <- qr.Q(qr(B))
+ B <- matrix(rnorm(N), nrow = m, ncol = n)
+ V <- qr.Q(qr(B))
+ d <- diag(sort(runif(50), decreasing = T)^6)
+ A <- U %*% d %*% t(V)
+ U1 <- svd(A)$u
+ d1 <- diag(svd(A)$d)
+ V1 <- svd(A)$v
+ normv <- function(x) max(svd(x)$d)
+ error.U1 <- normv(U1 - U)
+ error.V1 <- normv(V1 - V)
+ error.d1 <- normv(d1 - d)
+ error <- normv(U1 %*% d1 %*% t(V1) - A)
+ us <- which(abs(U1[1, ] - U[1, ]) > 10^-4)
+ vs <- which(abs(V1[1, ] - V[1, ]) > 10^-4)
+ j <- 1
+ for (j in seq_along(us)) U1[, us[j]] <- U1[, us[j]] * -1
+ for (j in seq_along(vs)) V1[, vs[j]] <- V1[, vs[j]] * -1
+ error.U1.f <- normv(U1 - U)
+ error.V1.f <- normv(V1 - V)
+ error.d1.f <- normv(d1 - d)
+ error.f <- normv(U1 %*% d1 %*% t(V1) - A)
+ res1 <- (c(error.U1, error.V1, error.d1, error))
+ res2 <- (c(error.U1.f, error.V1.f, error.d1.f, error.f))
+ res <- (rbind(res1, res2))
+ cond <- max(diag(d1))/min(diag(d1))
+ result <- list(res = res, cond = cond)
+ return(result)
+ }
> x <- svdtest1(123456)
> print(x$res)
[,1] [,2] [,3] [,4]
res1 2.000000e+00 2.000000e+00 4.440892e-16 2.917233e-15
res2 3.164195e-09 1.410472e-09 4.440892e-16 2.917233e-15
> print(x$cond)
[1] 87780167
> x <- svdtest1(12345)
> print(x$res)
[,1] [,2] [,3] [,4]
res1 2.00000e+00 2.000000e+00 7.771561e-16 1.098656e-15
res2 8.40177e-07 1.158272e-07 7.771561e-16 1.098656e-15
> print(x$cond)
[1] 234293764118
> x <- svdtest1(1234)
> print(x$res)
[,1] [,2] [,3] [,4]
res1 2.000000e+00 2.000000e+00 3.885781e-16 4.313957e-15
res2 7.787751e-08 5.416777e-08 3.885781e-16 4.313957e-15
> print(x$cond)
[1] 566641732936
> x <- svdtest1(123)
> print(x$res)
[,1] [,2] [,3] [,4]
res1 2.000000e+00 2.000000e+00 6.661338e-16 1.014207e-15
res2 5.975176e-08 4.430118e-08 6.661338e-16 1.014207e-15
> print(x$cond)
[1] 8279557898
> x <- svdtest1(12)
> print(x$res)
[,1] [,2] [,3] [,4]
res1 2.000000e+00 2.000000e+00 5.551115e-16 1.129351e-15
res2 3.210675e-09 3.756352e-09 5.551115e-16 1.129351e-15
> print(x$cond)
[1] 52212571023 |
Accuracy of individual U and V have drastically gone down.
This means that whichever algo you use, conditionality matters too Even if you use a robust svd algo, it does not mean anything.. If the matrix is illconditioned,then there are going to huge forward errors. One amazing aspect though is that the svd is remarkably backward stable.