3 Singular value decomposition

Singular value decomposition is a build-in function in R as well as P. Let us examine how long it takes to compute 1000 Hilbert matrices of dimension 9 and the corresponding SVD decomposition.

Here is the R code:

hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }

svd.test = function() {
X <- hilbert(9)[,1:6]
s <- svd(X)
D <- diag(s$d)
return(D)
}

svd.test.c = cmpfun(svd.test)

res <- benchmark( svd.test(), svd.test.c(),
                  columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"),
                  order="relative",
                  replications=1000)

The result is

         test replications elapsed relative user.self sys.self
2 svd.test.c()         1000   0.299 1.000000     0.283    0.016
1   svd.test()         1000   0.360 1.204013     0.284    0.014

R needs 0.36 seconds and the compiled version 0.299.

The P program looks a bit different, but not much:

source('stdlib.P')

hilbert.plus = function( x, y ) { var ; return( x+y ) }
hilbert = function( n ) {
  var i;
  i = 1:n
  return( 1/outer( i-1, i, hilbert.plus  ) )
}

svd.test = function() {
  var X, s, D;
  X <- hilbert(9)[,1:6]
  s <- svd(X)
  D <- diag(s$d)
  return(D)
}

repetitions = 1000
tm = time
i = 0
for ( i in (1:repetitions) ) val = svd.test()
print( time-tm )

Using P it takes only 0.106 seconds. This means, even in this, where most of time-consuming computations are done internally by the system, P runs 3.4 times faster than R and 2.82 times faster than R's byte compiler.

The Professional Edition needs 0.079 seconds. This is a speedup of 4.56 compared to R and 3.78 times faster than R's byte code.