[R] Positive Definite Matrix
Martin Maechler
maechler at stat.math.ethz.ch
Mon Jan 31 09:10:27 CET 2011
> 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
More information about the R-help
mailing list