[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