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.