On Jan 29, 2011, at 12:17 PM, Prof Brian Ripley wrote:
On Sat, 29 Jan 2011, David Winsemius wrote:
On Jan 29, 2011, at 10:11 AM, David Winsemius wrote:
On Jan 29, 2011, at 9:59 AM, John Fox wrote:
Dear David and Alex,
I'd be a little careful about testing exact equality as in all(M
== t(M) and
careful as well about a test such as all(eigen(M)$values > 0)
since real
arithmetic on a computer can't be counted on to be exact.
Which was why I pointed to that thread from 2005 and the existing
work that had been put into packages. If you want to substitute
all.equal for all, there might be fewer numerical false alarms,
but I would think there could be other potential problems that
might deserve warnings.
In addition to the two "is." functions cited earlier there is also
a "posdefify" function by Maechler in the sfsmisc package: "
Description : From a matrix m, construct a "close" positive
definite one."
But again, that is not usually what you want. There is no guarantee
that the result is positive-definite enough that the Cholesky
decomposition will work.
I don't see a Cholesky decomposition method being used in that
function. It appears to my reading to be following what would be
called an eigendecomposition.
--
David.
Give up on Cholesky factors unless you have a matrix you know must
be symmetric and strictly positive definite, and use the
eigendecomposition instead (setting negative eigenvalues to zero).
You can then work with the factorization to ensure that (for
example) variances are always non-negative because they are always
computed as sums of squares.
This sort of thing is done in many of the multivariate analysis
calculations in R (e.g. cmdscale) and in well-designed packages.
--
David.
On Jan 29, 2011, at 7:58 AM, David Winsemius wrote:
On Jan 29, 2011, at 7:22 AM, Alex Smith wrote:
Hello I am trying to determine wether a given matrix is
symmetric and
positive matrix. The matrix has real valued elements.
I have been reading about the cholesky method and another
method is
to find the eigenvalues. I cant understand how to implement
either of
the two. Can someone point me to the right direction. I have
used
?chol to see the help but if the matrix is not positive
definite it
comes up as error. I know how to the get the eigenvalues but
how can
I then put this into a program to check them as the just come
up with
$values.
Is checking that the eigenvalues are positive enough to
determine
wether the matrix is positive definite?
That is a fairly simple linear algebra fact that googling or
pulling
out a standard reference should have confirmed.
Just to be clear (since on the basis of some off-line
communications it
did not seem to be clear): A real, symmetric matrix is
Hermitian (and
therefore all of its eigenvalues are real). Further, it is
positive-
definite if and only if its eigenvalues are all positive.
qwe<-c(2,-1,0,-1,2,-1,0,1,2)
q<-matrix(qwe,nrow=3)
isPosDef <- function(M) { if ( all(M == t(M) ) ) { # first test
symmetric-ity
if ( all(eigen(M)$values > 0) )
{TRUE}
else {FALSE} } #
else {FALSE} # not symmetric
}
isPosDef(q)
[1] FALSE
m
[,1] [,2] [,3] [,4] [,5]
[1,] 1.0 0.0 0.5 -0.3 0.2
[2,] 0.0 1.0 0.1 0.0 0.0
[3,] 0.5 0.1 1.0 0.3 0.7
[4,] -0.3 0.0 0.3 1.0 0.4
[5,] 0.2 0.0 0.7 0.4 1.0
isPosDef(m)
[1] TRUE
You might want to look at prior postings by people more
knowledgeable than
me:
http://finzi.psych.upenn.edu/R/Rhelp02/archive/57794.html
Or look at what are probably better solutions in available
packages:
http://finzi.psych.upenn.edu/R/library/corpcor/html/rank.condition.html
http://finzi.psych.upenn.edu/R/library/matrixcalc/html/is.positive.definit
e.html
--
David.
this is the matrix that I know is positive definite.
eigen(m)
$values
[1] 2.0654025 1.3391291 1.0027378 0.3956079 0.1971228
$vectors
[,1] [,2] [,3] [,4] [,5]
[1,] -0.32843233 0.69840166 0.080549876 0.44379474
0.44824689
[2,] -0.06080335 0.03564769 -0.993062427 -0.01474690
0.09296096
[3,] -0.64780034 0.12089168 -0.027187620 0.08912912
-0.74636235
[4,] -0.31765040 -0.68827876 0.007856812 0.60775962
0.23651023
[5,] -0.60653780 -0.15040584 0.080856897 -0.65231358
0.42123526
and this are the eigenvalues and eigenvectors.
I thought of using
eigen(m,only.values=T)
$values
[1] 2.0654025 1.3391291 1.0027378 0.3956079 0.1971228
$vectors
NULL
m <- matrix(scan(textConnection("
1.0 0.0 0.5 -0.3 0.2
0.0 1.0 0.1 0.0 0.0
0.5 0.1 1.0 0.3 0.7
-0.3 0.0 0.3 1.0 0.4
0.2 0.0 0.7 0.4 1.0
")), 5, byrow=TRUE)
#Read 25 items
m
[,1] [,2] [,3] [,4] [,5]
[1,] 1.0 0.0 0.5 -0.3 0.2
[2,] 0.0 1.0 0.1 0.0 0.0
[3,] 0.5 0.1 1.0 0.3 0.7
[4,] -0.3 0.0 0.3 1.0 0.4
[5,] 0.2 0.0 0.7 0.4 1.0
all( eigen(m)$values >0 )
#[1] TRUE
Then i thought of using logical expression to determine if
there are
negative eigenvalues but couldnt work. I dont know what error
this is
b<-(a<0)
Error: (list) object cannot be coerced to type 'double'
??? where did "a" and "b" come from?
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Brian D. Ripley, [email protected]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
David Winsemius, MD
West Hartford, CT
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
