Hello thank you for so much input. I am afraid that I am fairly new to this
and some of the material is above my understanding for the time being and I
am not looking at anything very complex as you will see from the program I
will include.
I have been talking to a fellow colleague I am working with and have come
down to this.
Your input has been invaluable to me learning more about R thank you!
1.As an example, I do not wish to find a near pd matrix but I want to test
wether it is pd or not so I know if I want it or not,as you will see on the
program I want two outputs depending on if it is a pd and symmetric. So I
have written this for one symmetric matrix and it works and one for a non
symmetric one slightly changed but it doesnt work,if you could point out why
this is going wrong that would be helpful. It cant compare the complex
eigenvalues.
#symmetric matrix
elements<-c(1,0,0.5,-0.3,0.2,0,1,0.1,0,0,0.5,0.1,1,0.3,0.7,-0.3,0,0.3,1,0.4,0.2,0,0.7,3,1)
m<-matrix(elements,nrow=5,byrow=F)
sym<-isSymmetric(m)
eigen(m)
posdef1<-(all(eigen(m)$values > 0))
posdef1
m1<-m[c(1,2),c(1,2)]
m2<-m[c(1,2),c(3,4,5)]
m3<-m[c(3,4,5),c(1,2)]
m4<-m[c(3,4,5),c(3,4,5)]
im4<-solve(m4)
(m1-m2%*%im4%*%m3)*sym*posdef1
#
#not symmetric
elements<-c(1,0,0.5,-0.3,0.2,0,1,0.1,0,0,0.5,0.1,1,0.3,0.7,-0.3,0,0.3,1,0.4,0.2,0,0.7,3,1)
m<-matrix(elements,nrow=5,byrow=F)
sym<-isSymmetric(m)
sym
eigen(m)
posdef1<-(all(eigen(m)$values > 0))
posdef1
m1<-m[c(1,2),c(1,2)]
m2<-m[c(1,2),c(3,4,5)]
m3<-m[c(3,4,5),c(1,2)]
m4<-m[c(3,4,5),c(3,4,5)]
im4<-solve(m4)
(m1-m2%*%im4%*%m3)*sym*posdef1
#

2. My colleague said why dont we use the Sylvester's Criterion and test the
leading minors. Would something like that work?


On Mon, Jan 31, 2011 at 9:50 AM, Spencer Graves <
spencer.graves@structuremonitoring.com> wrote:

> 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@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@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@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
>
>

	[[alternative HTML version deleted]]

