Positive Definite Matrix
Hi, Martin: Thank you! (not only for your responses in this email thread but in helping create R generally and many of these functions in particular.) Spencer
On 1/31/2011 12:10 AM, Martin Maechler wrote:
> I think the bottom line can be summarized as
> follows:
> 1. Give up on Cholesky factors unless you have a
> matrix you know must be symmetric and strictly positive
> definite. (I seem to recall having had problems with chol
> even with matrices that were theoretically positive or
> nonnegative definite but were not because of round off.
> However, I can't produce an example right now, so I'm not
> sure of that.)
(other respondents to this thread mentioned such scenarios, they are not at all uncommon)
> 2. If you must test whether a matrix is
> summetric, try all.equal(A, t(A)). From the discussion, I
> had the impression that this might not always do what you
> want, but it should be better than all(A==t(A)). It is
> better still to decide from theory whether the matrix
> should be symmetric.
Hmm, yes: Exactly for this reason, R has had a *generic* function
isSymmetric()
-------------
for quite a while:
In "base R" it uses all.equal() {but with tightened default tolerance},
but e.g., in the Matrix package,
it "decides from theory" --- the Matrix package containing
quite a few Matrix classes that are symmetric "by definition".
So my recommendation really is to use isSymmetric().
> 3. Work with the Ae = eigen(A,
> symmetric=TRUE) or eigen((A+t(A))/2, symmetric=TRUE).
> From here, Ae$values<- pmax(Ae$values, 0) ensures that A
> will be positive semidefinite (aka nonnegative definite).
> If it must be positive definite, use Ae$values<-
> pmax(Ae$values, eps), with eps>0 chosen to make it as
> positive definite as you want.
hmm, almost: The above trick has been the origin and basic building block of posdefify() from the sfsmisc package, mentioned earlier in this thread. But I have mentioned that there are much better algorithms nowadays, and Matrix::nearPD() uses one of them .. actually with variations on the theme aka optional arguments.
> 4. To the maximum extent feasible, work with
> Ae, not A. Prof. Ripley noted, "You can then work with
> [this] 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."
yes, or---as mentioned by Prof Ripley as well---compute a
square root of the matrix {e.g. via the eigen() decomposition
with modified eigenvalues} and work with that.
Unfortunately, in quite a few situations you just need a
pos.def. matrix to be passed to another R function as
cov / cor matrix, and their, nearPD() comes very handy.
> Hope this helps. Spencer
It did, thank you, Martin
> On 1/30/2011 3:02 AM, Alex Smith wrote:
>> Thank you for all your input but I'm afraid I dont know
>> what the final conclusion is. I will have to check the
>> the eigenvalues if any are negative. Why would setting
>> them to zero make a difference? Sorry to drag this on.
>>
>> Thanks
>>
>> On Sat, Jan 29, 2011 at 9:00 PM, Prof Brian
>> Ripley<ripley at stats.ox.ac.uk>wrote:
>>
>>> On Sat, 29 Jan 2011, David Winsemius wrote:
>>>
>>>
>>>> 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.
>>>>
>>> Correct, but my point is that one does not usually want
>>> a
>>>
>>> "close" positive definite one
>>>
>>> but a 'square root'.
>>>
>>>
>>>
>>>> --
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?
-- Brian D. Ripley, ripley at stats.ox.ac.uk 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
-- Brian D. Ripley, ripley at stats.ox.ac.uk 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
Spencer Graves, PE, PhD President and Chief Operating Officer Structure Inspection and Monitoring, Inc. 751 Emerson Ct. San Jos?, CA 95126 ph: 408-655-4567