Skip to content

[Rcpp-devel] trans() changed in latest RcppArmadillo

9 messages · Conrad Sand, Savitsky, Terrance, Baptiste Auguie

Conrad Sand
#
On 30 May 2011, Terrance Savitsky wrote:
Hi, I'm the main author of Armadillo.

I'm interested in hearing about all regressions -- can you provide
more details ?

There have been a lot of changes between Armadillo 1.2 and the latest
beta (1.99.3).  I'd like to shake out all known bugs before releasing
2.0.


Cheers,
Conrad

--
Dr Conrad Sanderson : Sr Research Scientist : NICTA : http://arma.sf.net/cs
Savitsky, Terrance
# 
Hello Dr. Sanderson.  Thank you, so much, for your work; along with the
work of Dirk and Romain, your tools make it possible for me to have
impact I otherwise would never achieve.  Anyway, I separately tested my
(typical) use of trans() in a simple code and received no problems.  I
have various Rcpp-written functions that perform Bayesian regression
modeling; so I'm simulating samples from a set of conditional
distributions.   Until last night, I'd used these functions for some
months on a particular dataset without incident (though I've
intermittently tuned them for speed and functionality).  Anyway, further
analysis suggests that I'm experiencing numerical instability not
previously an issue.  I noted that Armadillo included some changes for
faster inverse computation on small size matrices.  Was there an
algorithm change that might allow such instability?  I intend to test by
replacing inv() with pinv().

Thanks, Terrance 

-----Original Message-----
From: Conrad Sand [mailto:conradsand.rcpp at gmail.com] 
Sent: Monday, May 30, 2011 6:34 AM
To: rcpp-devel at r-forge.wu-wien.ac.at
Cc: Savitsky, Terrance
Subject: Re: trans() changed in latest RcppArmadillo
On 30 May 2011, Terrance Savitsky wrote:
(0.2.20)
fix it.
my
Hi, I'm the main author of Armadillo.

I'm interested in hearing about all regressions -- can you provide
more details ?

There have been a lot of changes between Armadillo 1.2 and the latest
beta (1.99.3).  I'd like to shake out all known bugs before releasing
2.0.


Cheers,
Conrad

--
Dr Conrad Sanderson : Sr Research Scientist : NICTA :
http://arma.sf.net/cs

__________________________________________________________________________

This email message is for the sole use of the intended recipient(s) and
may contain confidential information. Any unauthorized review, use,
disclosure or distribution is prohibited. If you are not the intended
recipient, please contact the sender by reply email and destroy all copies
of the original message.
Savitsky, Terrance
# 
Dr. Sanderson, I've been able to verify that my issue resides in the
inv() function (that I typically apply to small p x p matrices, where p
= 3 - 10).  In particular, the new version/algorithm induces numerical
instability.  I've not yet tested that the results of the inv()
computation are generally accurate, only that they are not numerically
robust in comparison to 0.2.19.  I'm up against a deadline, so I
switched back to 0.2.19, which resolves my problem, for now, but will
provide a reproducible example when finished with my work.

Terrance  

-----Original Message-----
From: Conrad Sand [mailto:conradsand.rcpp at gmail.com] 
Sent: Monday, May 30, 2011 6:34 AM
To: rcpp-devel at r-forge.wu-wien.ac.at
Cc: Savitsky, Terrance
Subject: Re: trans() changed in latest RcppArmadillo
On 30 May 2011, Terrance Savitsky wrote:
(0.2.20)
fix it.
my
Hi, I'm the main author of Armadillo.

I'm interested in hearing about all regressions -- can you provide
more details ?

There have been a lot of changes between Armadillo 1.2 and the latest
beta (1.99.3).  I'd like to shake out all known bugs before releasing
2.0.


Cheers,
Conrad

--
Dr Conrad Sanderson : Sr Research Scientist : NICTA :
http://arma.sf.net/cs

__________________________________________________________________________

This email message is for the sole use of the intended recipient(s) and
may contain confidential information. Any unauthorized review, use,
disclosure or distribution is prohibited. If you are not the intended
recipient, please contact the sender by reply email and destroy all copies
of the original message.
Baptiste Auguie
# 
I was a bit optimistic yesterday; I have not yet been able to produce
a minimal example but it seems trans/strans was not the end of the
story in my code. Something else was broken with the change to
0.2.21.Thankfully the end results are very visibly wrong.

Best regards,

Baptiste
On 31 May 2011 07:02, Savitsky, Terrance <savitsky at rand.org> wrote:
Baptiste Auguie
# 
After a couple hours switching back and forth between versions, I
found a problem with the scalar product of complex vectors. Here is a
minimal example,

library(inline)
require( RcppArmadillo )
	
## let's calculate this product
c(-1-1i, 1-1i) %*% c(1i, -1i)
 ## 0-2i

## trans() with v0.2.19
fx <- cxxfunction( signature() , '
 arma::cx_colvec A = " (-1,-1)    (+1,-1) ;", B = "(+0,1)  (+0,-1) ;" ;
		arma::cx_colvec res = trans(A) * B;

		return wrap( res ) ;
	', plugin = "RcppArmadillo" )

fx()
## 0-2i

## strans() with v0.2.21
fx <- cxxfunction( signature() , '
 arma::cx_colvec A = " (-1,-1)    (+1,-1) ;", B = "(+0,1)  (+0,-1) ;" ;
		arma::cx_colvec res = strans(A) *B;

		return wrap( res ) ;
	', plugin = "RcppArmadillo" )

fx()
## 1-1i

Best regards,

baptiste
On 31 May 2011 10:01, baptiste auguie <baptiste.auguie at googlemail.com> wrote:
Baptiste Auguie
# 
Sorry, I had too many spaces in the vectors. The conclusion is the same though.

library(inline)
require( RcppArmadillo )
	
## let's calculate this product
c(-1-1i, 1-1i) %*% c(1i, -1i)
 ## 0-2i

## trans() with v0.2.19
fx <- cxxfunction( signature() , '
 arma::cx_colvec A = "(-1,-1) (+1,-1);", B = "(+0,1) (+0,-1);" ;
		arma::cx_colvec res = trans(A) * B;

		return wrap( res ) ;
	', plugin = "RcppArmadillo" )

fx()
## 0-2i

## strans() with v0.2.21
fx <- cxxfunction( signature() , '
 arma::cx_colvec A = "(-1,-1) (+1,-1);", B = "(+0,1) (+0,-1);" ;
		arma::cx_colvec res = strans(A) * B;
		return wrap( res ) ;
	', plugin = "RcppArmadillo" )

fx()
## 1-1i
On 31 May 2011 11:42, baptiste auguie <baptiste.auguie at googlemail.com> wrote:
Conrad Sand
#
On 31 May 2011 05:02, Savitsky, Terrance <savitsky at rand.org> wrote:
This is a bit strange, as the underlying algorithm for inv() hasn't
changed.  For matrix sizes <= 4x4, a fast version is used (and has
been since the early days of Armadillo).  For sizes >= 5x5, Lapack is
used.

Can you provide an example matrix which provides different results
with the new Armadillo ?

btw, in the new Armadillo the solve() function now calls inv() for
matrix sizes <= 4x4.  Otherwise it uses Lapack.
That would be great -- as the issue affects you, it's likely to affect
other people.


Cheers,
Conrad

--
Dr Conrad Sanderson : Sr Research Scientist : NICTA : http://arma.sf.net/cs
Conrad Sand
# 
Hi Baptiste,

Thanks for the bug report.  I'll take a look at the underlying issues.
On 31 May 2011 09:48, baptiste auguie <baptiste.auguie at googlemail.com> wrote:
Conrad Sand
# 
I've found the cause of the issue (involves handling of small matrices).

The fix is easy (and already done in the SVN repo), but I currently
don't have the time to roll out another release.  I'll aim to release
a fix on the weekend, which will also give me time to double-check if
there is a problem in inv().

In the meantime I recommend sticking to the previous version of
Armadillo and RcppArmadillo.

btw, for dot products I recommend using the dot() and cdot() functions
-- they're generally faster than going through the multiplication
operator.
On 31 May 2011 11:39, Conrad Sand <conradsand.rcpp at gmail.com> wrote: