You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
</style></head><body><divclass="content"><h2>Contents</h2><div><ul><li><ahref="#2">nearestSPD works on any matrix, and it is reasonably fast.</a></li><li><ahref="#7">A realistic test case</a></li></ul></div><preclass="codeinput"><spanclass="comment">% Neearest Symmetric, Positive Definite matrices</span>
69
+
<spanclass="comment">%</span>
70
+
<spanclass="comment">% This tool saves your covariance matrices, turning them into something</span>
71
+
<spanclass="comment">% that really does have the property you will need. That is, when you are</span>
72
+
<spanclass="comment">% trying to use a covariance matrix in a tool like mvnrnd, it makes no</span>
73
+
<spanclass="comment">% sense if your matrix is not positive definite. So mvnrnd will fail in</span>
74
+
<spanclass="comment">% that case.</span>
75
+
<spanclass="comment">%</span>
76
+
<spanclass="comment">% But sometimes, it appears that users end up with matrices that are NOT</span>
77
+
<spanclass="comment">% symmetric and positive definite (commonly abbreviated as SPD) and they</span>
78
+
<spanclass="comment">% still wish to use them to generate random numbers, often in a tool like</span>
79
+
<spanclass="comment">% mvnrnd. A solution is to find the NEAREST matrix that has the desired</span>
80
+
<spanclass="comment">% property of being SPD.</span>
81
+
<spanclass="comment">%</span>
82
+
<spanclass="comment">% I see the question come up every once in a while, so I looked in the file</span>
83
+
<spanclass="comment">% exchange to see what is in there. All I found was nearest_posdef. While</span>
84
+
<spanclass="comment">% this usually almost works, it could be better. It actually failed</span>
85
+
<spanclass="comment">% completely on most of my test cases, and it was not as fast as I would</span>
86
+
<spanclass="comment">% like, using an optimization. In fact, in the comments to nearest_posdef,</span>
87
+
<spanclass="comment">% a logical alternative was posed. That alternative too has its failures,</span>
88
+
<spanclass="comment">% so I wrote nearestSPD.</span>
89
+
</pre><h2>nearestSPD works on any matrix, and it is reasonably fast.<aname="2"></a></h2><p>As a test, randn generates a matrix that is not symmetric nor is it at all positive definite in general.</p><preclass="codeinput">U = randn(100);
90
+
</pre><p>nearestSPD will be able to convert U into something that is indeed SPD, and for a 100 by 100 matrix, do it quickly enough</p><preclass="codeinput">tic,Uj = nearestSPD(U);toc
91
+
</pre><preclass="codeoutput">Elapsed time is 0.005662 seconds.
92
+
</pre><p>The ultimate test of course, is to use chol. If chol returns a second argument that is zero, then MATLAB (and mvnrnd) will be happy!</p><preclass="codeinput">[R,p] = chol(Uj);
93
+
p
94
+
</pre><preclass="codeoutput">p =
95
+
0
96
+
</pre><p>As you can see, mvnrnd did not complain at all.</p><preclass="codeinput">mvnrnd(zeros(1,100),Uj,1)
</pre><p>nearest_posdef would have failed here, as U was not even symmetric, nor does it even have positive diagonal entries.</p><h2>A realistic test case<aname="7"></a></h2><p>Next I'll try a simpler test case. This one will have positive diamgonal entries, and it will indeed be symmetric. So this matrix is much closer to a true covariance matrix than that first mess we tried. And since nearest_posdef was quite slow on a 100x100 matrix, I'll use something smaller.</p><preclass="codeinput">U = rand(25);
129
+
U = (U + U')/2;
130
+
</pre><p>Really, it is meaningless as a covariance matrix, because it is clearly not positive definite. We can see many negative eigenvalues, and chol gets upset. So mvnrnd would fail here.</p><preclass="codeinput">eig(U)'
</pre><p>nearest_posdef took a bit of time, about 9 seconds on my machine. Admittedly, much of that time was wasted in doing fancy graphics that nobody actually needs if they just need a result.</p><preclass="codeinput">tic,Um = nearest_posdef(U);toc
145
+
</pre><preclass="codeoutput">Elapsed time is 8.768167 seconds.
146
+
</pre><imgvspace="5" hspace="5" src="nearestSPD_demo_01.png" alt=""><p>Is Um truly positive definite according to chol? Sadly, it is usually not.</p><preclass="codeinput">[R,p] = chol(Um);
147
+
p
148
+
</pre><preclass="codeoutput">p =
149
+
18
150
+
</pre><p>We can see how it failed, by looking at what eig returns.</p><preclass="codeinput">eig(Um)
151
+
</pre><preclass="codeoutput">ans =
152
+
11.8426 + 0.0000i
153
+
1.7845 + 0.0000i
154
+
1.5092 + 0.0000i
155
+
1.2512 + 0.0000i
156
+
1.0957 + 0.0000i
157
+
0.7042 + 0.0000i
158
+
0.4589 + 0.0000i
159
+
0.3627 + 0.0000i
160
+
0.3080 + 0.0000i
161
+
0.2593 + 0.0000i
162
+
0.1571 + 0.0000i
163
+
0.0417 + 0.0000i
164
+
0.0387 + 0.0000i
165
+
0.0290 + 0.0000i
166
+
0.0190 + 0.0000i
167
+
0.0052 + 0.0000i
168
+
0.0021 + 0.0000i
169
+
-0.0000 + 0.0000i
170
+
-0.0000 - 0.0000i
171
+
-0.0000 + 0.0000i
172
+
0.0000 + 0.0000i
173
+
0.0000 + 0.0000i
174
+
-0.0000 + 0.0000i
175
+
-0.0000 - 0.0000i
176
+
0.0000 + 0.0000i
177
+
</pre><p>There will usually be some tiny negative eigenvalues</p><preclass="codeinput">min(real(eig(Um)))
178
+
</pre><preclass="codeoutput">ans =
179
+
-1.8002e-16
180
+
</pre><p>and sometimes even some imaginary eigenvalues. All usually tiny, but still enough to upset chol.</p><preclass="codeinput">max(imag(eig(Um)))
181
+
</pre><preclass="codeoutput">ans =
182
+
2.2548e-16
183
+
</pre><p>The trick suggested by Shuo Han is pretty fast, but it too fails. Since U is already symmetric, we need not symmetrize it first, but chol still gets upset.</p><p>Note that the slash used by Shuo was not a good idea. transpose would have been sufficient.</p><preclass="codeinput">[V,D] = eig(U);
184
+
U_psd = V * max(D,0) / V;
185
+
[R,p] = chol(U_psd);
186
+
p
187
+
</pre><preclass="codeoutput">p =
188
+
13
189
+
</pre><p>Whereas nearestSPD works nicely.</p><preclass="codeinput">Uj = nearestSPD(U);
190
+
[R,p] = chol(Uj);
191
+
</pre><p>nearestSPD returns a solution that is a bit closer to the original matrix U too. Thus comparing 1,2,inf and Frobenious norms, nearestSPD was better under all norms in my tests, even though it is designed only to optimize the Frobenious norm.</p><preclass="codeinput">[norm(U - Um,1), norm(U - Um,2), norm(U - Um,inf), norm(U - Um,<spanclass="string">'fro'</span>)]
0 commit comments