[R] Positive Definite Matrix
Spencer Graves
spencer.graves at structuremonitoring.com
Mon Jan 31 10:50:12 CET 2011
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
More information about the R-help
mailing list