s-news
[Top] [All Lists]

Re: chol error in Finmetric's SUR

To: Eric Gravengaard <gravengaard@gmail.com>, "Stahel, Christof" <stahel_2@cob.osu.edu>
Subject: Re: chol error in Finmetric's SUR
From: Tony Plate <tplate@blackmesacapital.com>
Date: Mon, 19 Apr 2004 15:39:42 -0600
Cc: s-news@lists.biostat.wustl.edu
In-reply-to: <EFD4B5D2.1CEC9D8@mail.gmail.com>
References: <1C0CCEE45B00E74D8FBBB1AE6A472E855D66E3@COBMAIL.ad.cob.ohio-state.edu> <EFD4B5D2.1CEC9D8@mail.gmail.com>
The undocumented S-plus function choleski() does not have the problem that chol() has with matrices containing minor deviations from symmetry. Although it returns data in a different format than chol(), it's pretty easy to use the results, especially if no pivoting was performed.

I first discovered this when using a new Windows Xeon workstation. The fast BLAS libraries can return a non-symmetric matrix as the result of crossprod(x) (== t(x)%*%x). If this is the source of your problem, doing intelmkl.use(set=F) (on an Intel Windows system) may help, though it will probably increase run times for matrix operations.

The following old message from Insightful support may help (however, note that there might be problems with their suggested solution if the result of choleski() does not have full rank -- I think in that case the result of choleski() may need pivoting to be the same as the corresponding result from chol()).

hope this helps,

Tony Plate

Subject: RE: t(x)%*%x is non-symmetric in S-plus 6.1 on HP xw8000
Date: Thu, 27 Mar 2003 08:08:50 -0800

Dear Tony Plate,

Thank you for your question regarding choleski().  I created the
enhancement request to have an online help file created for
this function, with links to the chol() function.

You can use the following to convert the output of choleski
into something like the output of chol().  solve.upper wants it's input
to be a matrix with an attribute called "rank", while choleski's
ouput is a list, with an element for the matrix and an element for
the rank (and some other things).

f <- function(x) {
        # expect x to be output of choleski(), with names
        # "R"         "rank"      "condition" "pivot"
        # solve.upper wants a matrix with an extra attribute "rank"
        val <- x$R
        attr(val, "rank") <- x$rank
        oldClass(val) <- "upper"
        val
}

I hope this helps.  Please let me know if you have additional questions.

Best Regards,

Sherry LaMonica
Insightful Technical Support


At Monday 02:21 PM 4/19/2004, Eric Gravengaard wrote:
Christof,

I've had similar errors when using a Choleski decomposition. It comes
down to the matrix not being symmetric. The problems usually arise
when the matrix is symmetric, but because of machine precision
limitations doesn't look symmetric to a machine. I do not have a work
around except to use Matlab to do this because their implementation
does not have this limitation.

Do others have work arounds or possibly another implementation of the
Choleski decomposition that would help to solve this problem?

-Eric

On Sat, 17 Apr 2004 12:49:12 -0400, Stahel, Christof
<stahel_2@cob.osu.edu> wrote:
>
> Dear all,
> I receive the following error
>
> Problem in chol(xpx[X.idx[[i]], X.idx[[i]]]): Matrix must be symmetric
> Use traceback() to see the call stack
> when I call SUR in Finmetrics using the sample code straight from the help file. Any idea what causes the problem?
>
> # simulate a multivariate data set.
> set.seed(10)
> x1 <- abs(rnorm(100))
> x2 <- x1^2
> sim.dat <- data.frame(x1=x1,
>                      x2=x2,
>                      y1=1+2*x1+rnorm(100),
>                      y2=4+5*x2+rnorm(100))
> # estimate an SUR model
> SUR(list(y1~x1, y2~x2), data=sim.dat)
>
> I am using  release 6.0.2 for windows on a 3.06 Pentium 4 HT running XP.
>
> Thanks VERY much for your help & input!
> C
> --------------------------------------------------------------------
> This message was distributed by s-news@lists.biostat.wustl.edu.  To
> unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
> the BODY of the message:  unsubscribe s-news
>
--------------------------------------------------------------------
This message was distributed by s-news@lists.biostat.wustl.edu.  To
unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
the BODY of the message:  unsubscribe s-news


<Prev in Thread] Current Thread [Next in Thread>