1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043

<!DOCTYPE HTML>
<HEAD><TITLE>Surface Evolver Documentation  Stability Tutorial</title>
<link rel="stylesheet" type="text/css" href="evdocstyle.css" />
<meta httpequiv="ContentType" content="text/html;charset=utf8" >
</head>
<BODY>
<!NewPage>
<h1 class="center">
<a href="http://www.susqu.edu/brakke/evolver/evolver.htm" class="comic">
Surface Evolver</a> Documentation</h1>
<a href="evolver.htm#doctop">Back to top of Surface Evolver documentation.</a>
<a href="index.htm">Index.</a>
<a id="hessiantutorial"></a>
<h1>Hessian, Eigenvalues, Eigenvectors, and Stability
in the Surface Evolver</h1>
Our motto: The purpose of calculus is to turn every problem
into a linear algebra problem.
<p>
Contents:
<ol>
<li> <a href="#configurationspace">How to think of surfaces discretized in the Evolver.</a>
<li> <a href="#hessianintro">What is the Hessian?</a>
<li> <a href="#eigenvalues">Eigenvalues and eigenvectors.</a>
<li> <a href="#eigenpairinterpretation">Physical interpretation of eigenvalues and eigenvectors.</a>
<li> <a href="#eigenprobe">Eigenprobe command.</a>
<li> <a href="#ritz">Ritz command.</a>
<li> <a href="#square">The Square example.</a>
<li> <a href="#hessianmetric">Hessian with metric.</a>
<li> <a href="#hessiannormalmode">Normal motion mode.</a>
<li> <a href="#eigenvectorvisualization">Visualizing eigenvectors.</a>
<li> <a href="#hessianiteration">Hessian iteration.</a>
<li> <a href="#hessian_seek">Hessian_seek.</a>
<li> <a href="#saddleexample">Saddle.</a>
<li> <a href="#unstablesurfaces">Unstable surfaces.</a>
<li> <a href="#eigenvalueaccuracy">Accuracy of eigenvalues.</a>
<li> <a href="#hessiancaveats">Caveats and warnings.</a>
</ol>
For background,
see any decent linear algebra text, such as Strang <a href="biblio.htm#refSG">[SG]</a>, especially chapter 6.
For more on stability and bifurcations, see Arnol'd <a href="biblio.htm#refAV">[AV]</a> or Hale and
Kocak <a href="biblio.htm#refHK">[HK]</a>.
<hr><a id="configurationspace"></a><h2>
1. How to think of surfaces discretized in the Evolver.</h2>
The Surface Evolver parameterizes a surface in terms of vertex
coordinates. Restrict attention to one fixed triangulation
of a surface. Let X be the vector of all vertex coordinates,
containing N entries, where N could be thousands. Call X the
"configuration vector", and call the Ndimensional vector space
it inhabits "configuration space". The total energy E(X) is
a scalar function on configuration space. Visualize the graph
of the energy function as a rolling landscape with mountains, valleys,
and mountain passes. The gradient of the energy is the steepest
uphill direction. The ordinary
'<a href="single.htm#g" class="keyword">g</a>' iteration step of Evolver
minimizes energy by moving in the downhill direction, the negative
of the gradient direction, seeking a local minimum.
<p>
A <em>constraint</em> is a scalar function on configuration space
such that X is restricted to lie on a particular level set. A
global constraint, such
as a volume constraint, determines a hypersurface, that is, it removes one
degree of freedom from the configuration space.
<a href="constrnt.htm#levelsetconstraints">Levelset constraints</a>
on vertices,
such as x^2 + y^2 = 4, are actually a separate constraint for each vertex
they apply to. Each removes one degree of freedom for each vertex
it applies to. The intersection of all the constraint hypersurfaces
is called the "feasible set", and is the set of all possible
configurations satisfying all the constraints. If there are any nonlinear
constraints, such as curved walls or volume constraints, then the
feasible set will be curved. When a surface is moved, by
'<a href="single.htm#g" class="keyword">g</a>' for
example, the vector for the direction of motion is always first
projected to be tangent to the feasible set. However, if the
feasible set is curved, straightline motion will deviate
slightly from the feasible set. Therefore, after every motion
X is projected back to the feasible set using Newton's Method.
<hr><a id="hessianintro"></a><h2>
2. What is the Hessian?</h2>
Consider a particular surface in a particular configuration X,
and consider a small perturbation Y added to X. Then the energy
may be expanded in a Taylor series:
<pre>
E(X+Y) = E<sub>0</sub> + G<sup>T</sup> Y + 1/2 Y<sup>T</sup> H Y + ...
</pre>
Here the constant E<sub>0</sub> is the original energy, E<sub>0</sub> = E(X).
G is the vector of first derivatives, the gradient, which is a
vector the same dimension as Y or X (G<sup>T</sup> is a 1form, if you want
to get technical).
Here <sup>T</sup> denotes the transpose.
H is the square matrix of second
derivatives, the Hessian.
The Hessian
is automatically a symmetric matrix. The gradient G determines the
best linear approximation to the energy, and the Hessian H determines
the best quadratic approximation.
<p>
<b>Positive definiteness.</b>
If the quadratic term 1/2 Y<sup>T</sup> H Y is always positive for
any nonzero Y, then H is said to be <em> positive definite</em>.
At an equilibrium point, this means the point is a strict local minimum;
this is the multivariable version of the second derivative test for
a minimum. If the quadratic term is zero or positive, then
H is <em>positive semidefinite</em>. No conclusion can be drawn
about being a local minimum, since thirdorder terms may prevent that.
If the quadratic term has positive and negative values, then H is
<em>indefinite</em>.
<p>
<b>Constraints.</b> If there are constraints on the surface,
then the the perturbation
vector Y is required to be tangent to the feasible set in configuration
space. The curvature of the feasible set can have an effect on H,
but this is taken care of automatically by the Evolver. This effect
arises from the slight change in energy due to the projection back
to the feasible set after motion. This does not affect the gradient,
but it does affect the Hessian. In particular, if Q is the Hessian
of a global constraint and q is the Lagrange multiplier for the constraint,
then the adjusted energy Hessian is H  qQ. Therefore it is necessary
to do a '<a href="single.htm#g" class="keyword">g</a>' step to calculate Lagrange
multipliers before doing any Hessians when there are global constraints.
<hr><a id="eigenvectors"></a><h2>
<a id="eigenvalues"></a>
3. Eigenvalues and eigenvectors.</h2>
The Hessian H is a real symmetric matrix. Therefore it can
be diagonalized by an orthogonal change of basis of configuration space.
The new basis vectors are called <em>eigenvectors</em>, and the entries
on the diagonal version of H are called <em>eigenvalues</em>.
In this eigenvector
basis, the shape of the graph of the quadratic term becomes obvious.
Directions
along eigenvectors with negative eigenvalues have curvature
downward. Directions with positive eigenvalues have upward
curvature. Remember that eigenvectors are in configuration
space, so each eigenvector represents a particular perturbation
of the surface.
<p>
Another characterization of an eigenvector Q with an eigenvalue q
is
<pre> H Q = q Q.
</pre>
which is obvious in an eigenvector basis, but serves to define
eigenvectors in an arbitrary basis.
<p>
<a id="index"></a>
<b>Index and inertia.</b>
The number of negative eigenvalues of H is called the <em>index</em> of H.
<a id="inertia"></a>
The triple of numbers of negative, zero, and positive eigenvalues is called
the <em>inertia</em> of H. One can also speak of the inertia relative
to a value c as the inertia of HcI, i.e. the number of eigenvalues
less than c, equal to c, and greater than c.
Sylvester's Law of Inertia says that if P is a nonsingular matrix
of the same dimension as H, then the matrix
<pre>
K = P H P<sup>T</sup>
</pre>
has the same inertia as H, although not the same eigenvalues. The
Evolver can factor a Hessian into a form
<pre>
H = L D L<sup>T</sup>
</pre>
where L is a lower triangular matrix and D is diagonal. Hence the
inertia of H can be found by counting signs on the diagonal of D.
Inertia relative to c can be found by factoring HcI.
<p>
The eigenvectors associated to a particular eigenvalue form a
subspace, whose dimension is called the "multiplicity" of the
eigenvalue.
<b> With volume constraints.</b> If a surface is subject to
volume constraints or named quantity constraints, these are taken
into account in calculating eigenvalues and eigenvectors. In particular,
when defining a perturbation by a vectorfield, projection back to
the constraints is done at the end of the projection. So for a
nonlinear volume constraint (which is usually the case), adding
a volume constraint is different than taking the eigenvalues of
the unconstrained hessian relative to a linear subspace. It can
happen that adding a constraint can have a lower lowest eigenvalue
than without the constraint.
<hr><a id="eigenpairinterpretation"></a><h2>
4. Physical interpretation of eigenvalues and eigenvectors.</h2>
One can think of an eigenvector as a mode of pertubation of the surface,
and the corresponding eigenvalue an indicator of its stability.
Consider a surface at equilibrium. Recall that the gradient after a
small perturbation from equilibrium is
<pre>
grad E = H Y.
</pre>
So if we move along an eigenvector direction,
<pre>
grad E = H Q = q Q.
</pre>
Note the force is aligned exactly with the perturbation. Since
force is negative gradient, if q is positive the force will
restore the surface to equilibrium (stable perturbation),
and if q is negative the
force will cause the pertubation to grow (unstable perturbation).
Each eigenvector thus
represents a mode of perturbation. Furthermore, different modes
grow or decay independently, since in an eigenvector basis the
Hessian is diagonal and there are no terms linking the different
modes.
<hr><a id="eigenprobe"></a>
<h2>5. Eigenprobe command</h2>
The inertia of the Hessian with respect to some value
may be found with the <a href="commands.htm#eigenprobe" class="keyword">eigenprobe</a>
command.
The syntax to find the inertia relative to a probe value c is
"<code class="shortcode">eigenprobe c</code>".
For example, "<code class="shortcode">eigenprobe 1</code>"
would report the number of eigenvalues less than 1, equal to 1,
and greater than 1. For example, suppose we use cube.fe, and
do "<code class="shortcode">r; r; g 5; eigenprobe 1</code>". Then we get
<pre>
Vertices: 50 Edges: 144 Facets: 96 Bodies: 1 Facetedges: 288
Element memory: 39880
Total data memory: 294374 bytes.
Enter command: r
Vertices: 194 Edges: 576 Facets: 384 Bodies: 1 Facetedges: 1152
Element memory: 157384
Total data memory: 422022 bytes.
Enter command: g 5
5. area: 5.35288108682446 energy: 5.35288108682446 scale: 0.233721
4. area: 5.14937806689761 energy: 5.14937806689761 scale: 0.212928
3. area: 5.04077073041990 energy: 5.04077073041990 scale: 0.197501
2. area: 4.97600054126748 energy: 4.97600054126748 scale: 0.213378
1. area: 4.93581733363156 energy: 4.93581733363156 scale: 0.198322
Enter command: eigenprobe 1
Eigencounts: 18 <, 0 ==, 175 >
</pre>
This means the Hessian has 18 eigenvalues less than 1, none equal to 1,
and 175 greater than 1. Note the total number of eigenvalues is equal
to the total number of degrees of freedom of the surface, in this
case
<pre>
(# vertices)(# constraints) = 194  1 = 193 = 18+175
</pre>
where the constraint is the volume constraint. Why there is one
degree of freedom per vertex instead of three is explained below in
the <a href="#hessiannormalmode" class="keyword">normal_motion</a> section.
<p>
<hr>
<a id="ritz"></a>
<h2>6. Ritz command</h2>
For calculation of eigenvalues near a particular probe value,
there is the <a href="commands.htm#ritz" class="keyword">ritz</a>
command (named after the RayleighRitz algorithm). The syntax is
"<code>ritz(<em>c,n</em>)</code>", where c is the probe value and n is the desired
number of eigenvalues.
<p>
For example, evolve cube.fe with "<code class="shortcode">g 5; r; g 5; r; g 5</code>",
and do "<code class="shortcode">ritz(0,5)</code>".
You should get output like this:
<pre>
Enter command: ritz(0,5)
Eigencounts: 0 <, 0 ==, 193 >
1. 0.001687468263013
2. 0.001687468263013
3. 0.001687468263013
4. 0.2517282974725
5. 0.2517282974731
Iterations: 710. Total eigenvalue changes in last iteration: 4.0278891e013
</pre>
The first line of output is the inertia, exactly as in eigenprobe. Then
come the eigenvalues as the algorithm finds them, usually in order away
from the probe value.
<p>
How ritz works: For <code class="shortcode">ritz(c,n)</code>, Evolver creates a random
ndimensional subspace and applies shifted inverse iteration to it,
i.e. repeatedly applies (H  cI)<sup>1</sup> to the subspace and
orthonormalizes it. Eigenvalues of H near p correspond to large
eigenvalues of (H  cI)<sup>1</sup>, so the corresponding eigenvector
components grow by large factors each inverse iteration.
Evolver reports eigenvalues as they converge to machine precision. When all
desired eigenvalues converge (as by the total change in eigenvalues in an
iteration being less than 1e10), the iteration ends and the values are
printed. Lesser precision is shown on these to indicate they are not
converged to machine precision. You can interrupt ritz by hitting the
interrupt key (usually CTRLC), and it will report the current eigenvalue
approximations.
<p>NOTE: It is best to use probe values below the
least eigenvalue, so HcI is positive
definite. Makes the factoring more numerically stable.
<hr> <a id="square"></a>
<h2>7. The Square example</h2>
<img src="square.gif" alt="square"><br>
We now explore the eigenvalues of a very simple surface, a flat square
with fixed sides. For a continuous surface, classical calculus of variations
reveals that the Hessian for area is the Laplace operator. For
a square of side length pi, this has eigenvalues m<sup>2</sup> + n<sup>2</sup>
for m,n = 1,2,3,..., or in ascending sequence
2, 5, 5, 8, 10, 10, 13, 13, 17, 17, ..
(you may remember this from PDE analysis of the heat
equation). Run the datafile square.fe and refine it twice (no need to
evolve, since it is already flat). Do "<code class="shortcode">ritz(0,10)</code>".
<pre>
Enter command: ritz(0,10)
Eigencounts: 0 <, 0 ==, 25 >
1. 0.585786437626905
2. 1.386874070247247
3. 1.386874070247247
4. 2.000000000000000
5. 2.585786437626905
6. 2.585786437626906
7. 2.917607799707606
8. 2.9176077997076
9. 3.4142135623733
10. 4.0000000000000
</pre>
Doesn't look exactly like we expected. It does have multiplicity two
eigenvalues in the right places, though. Maybe we haven't refined enough.
Refine again and do "<code class="shortcode">ritz(0,10)</code>". You get
<pre>
Enter command: ritz(0,10)
Eigencounts: 0 <, 0 ==, 113 >
1. 0.152240934977427
2. 0.375490214588449
3. 0.375490214588449
4. 0.585786437626905
5. 0.738027372604331
6. 0.738027372604331
7. 0.927288973154334
8. 0.927288973154335
9. 1.225920309355705
10. 1.2259203093583
</pre>
Things are getting worse! The eigenvalues all shrunk drastically! They
are not converging! What's going on? The next section explains.
<hr><a id="hessianmetric"></a><h2>
8. Hessian with metric.</h2>
One would expect that refining the surface would lead the eigenvalues
to converge to the eigenvalues for the smooth surface. But as the
previous section showed, refining caused the eigenvectors to all
shrink by about a factor of 4. This is no
way to converge. The explanation is that so far we have been
looking at eigenvectors in slightly the wrong way. An eigenvector
is supposed to represent a mode of perturbation that is proportional
to the force. However, the response of a surface to a force need
not be numerically equal to the force. After all, forces and
displacements are different kinds of things. They have different
units: displacement has units of distance, and force has units
of energy per distance. In other words, displacement is a vector
and force is a covector. Note that matrix multiplication by
the Hessian H turns a vector into
a covector. In general, to turn a vector into an equivalent covector, one
needs an inner product, or metric. So far we have been using the Euclidean
inner product on configuration space, but that is not the proper
one to use if you want to approximate smooth surfaces. The
proper inner product of perturbations f and g of a smooth surface
is the integral over the surface of the pointwise inner product:
<pre>
/
<f,g> =  <f(x),g(x)> dA.
/
</pre>
In discrete terms, the inner product takes the form of a square
matrix M such that <Y,Z> = Y<sup>T</sup> M Z for vectors Y,Z. We want
this inner product to approximate integration with respect to
area. Having such an M, the eigenvalue equation becomes
<pre>
H Q = q M Q.
</pre>
Officially, Q is now called a "generalized eigenvector" and q is a
"generalized eigenvalue". But we will drop the "generalized".
An intuitive interpretation of this equation is as Newton's Law
of Motion,
<pre> Force = Mass * Acceleration </pre>
where HQ is the force, M is the mass, and qQ is the acceleration.
In other words, in an eigenmode, the acceleration is proportional
to the perturbation.
<p>
The Evolver toggle command
"<a href="toggle.htm#linear_metric" class="keyword">linear_metric</a>"
includes M in the eigenvector
calculations. Two metrics are available. In the simplest,
the "star metric", M is a diagonal matrix, and the "mass" associated with
each vertex is 1/3 of the area of the adjacent facets (1/3 since each
facet gets shared among 3 vertices). The other, the "linear metric",
extends functions from vertices to facets by linear interpolation
and integrates with respect to area. By default, "<code>linear_metric</code>"
uses a 50/50 mix of the two, which seems to work best. See the
more detailed discussion below in the <a href="#eigenvalueaccuracy">
eigenvalue accuracy</a> section. The fraction of linear metric
can be set by assigning the fraction to the internal variable
<a href="syntax.htm#linear_metric_mix" class="keyword">linear_metric_mix</a>.
By default, <code>linear_metric_mix</code> is 0.5. In quadratic mode, however,
only the quadratic interpolation metric is used, since the star metric
restricts convergence to order h<sup>2</sup> while the quadratic interpolation
metric permits h<sup>4</sup> convergence.
<p>
Example: Run square.fe, and refine twice. Do "<code class="code">linear_metric</code>" and
"<code class="code">ritz(0,10)</code>".
<pre>
Enter command: linear_metric
Using linear interpolation metric with Hessian.
Enter command: ritz(0,10)
Eigencounts: 0 <, 0 ==, 25 >
1. 2.036549240354335
2. 4.959720306550875
3. 4.959720306550875
4. 7.764634143334637
5. 10.098798069316683
6. 10.499717069102575
7. 12.274525789880887
8. 12.274525789880890
9. 15.7679634642721
10. 16.7214142904405
Iterations: 127. Total eigenvalue changes in last iteration: 1.9602375e014
</pre>
After refining again:
<pre>
Enter command: ritz(0,10)
Eigencounts: 0 <, 0 ==, 113 >
1. 2.009974216370147
2. 4.999499042685446
3. 4.999499042685451
4. 7.985384943789077
5. 10.090156419079085
6. 10.186524915471155
7. 13.008227978344527
8. 13.008227978344527
9. 17.242998229925817
10. 17.2429982299600
Iterations: 163. Total eigenvalue changes in last iteration: 1.9186042e014
</pre>
This looks much more convergent.
<p>
Using <code>linear_metric</code> does NOT change the inertia of the Hessian, by
Sylvester's Law of Inertia. So the moral of the story is that if
you care only about stability, you don't need to use <code>linear_metric</code>.
If you do care about the actual values of eigenvectors, and want them
relatively independent of your triangulation, then use <code>linear_metric</code>.
<hr><a id="hessiannormalmode"></a><h2>
9. Normal motion mode. </h2>
The alert reader will have notice that in the examples so far there
has been only one degree of freedom per vertex, instead of the three
one might expect, since a vertex has three degrees of freedom to
move in space. The answer is that in Hessian activities, it is usually
best to only allow motion perpendicular to the surface, and suppress
the two degrees of freedom of motion of a vertex tangential to the
surface. The reason is that tangential motion changes the energy
of the surface very little if at all, leading to many small eigenvalues
and a severely singular Hessian. For example, run square.fe, refine twice,
and do "<code class="shortcode">hessian_normal off</code>" to
enable all degrees of freedom. Now
"<code class="shortcode">eigenprobe 0</code>" reveals 50 zero eigenvalues:
<pre>Enter command: hessian_normal off
hessian_normal OFF. (was on)
Enter command: eigenprobe 0
Eigencounts: 0 <, 50 ==, 25 >
</pre>
For a curved surface, the extra eigenvalues generally won't be zero, but
they will all be close to zero, both positive and negative, which can
really foul things up. The default value of <code>hessian_normal</code> is therefore
the ON state. The moral of the story is to always leave hessian normal
on, unless you really really know what you are doing.
<p>
On some surfaces, for example soap films with triple junctions and
tetrahedral points, there are vertices with no natural normal vector.
Evolver <code>hessian_normal</code> mode assigns those points normal subspaces
instead, so that vertices on a triple line can move in a twodimensional
normal space perpendicular to the triple line, and tetrahedral points
can move in all three dimensions.
<p>
The reason for possible negative extra eigenvalues when hessian_normal
is off is that one rarely has the best possible vertex locations
for a given triangulation of
a surface, even when its overall shape is very close to optimal.
Vertices always seem to want to slither sideways to save very very
small amounts of energy. The amount of energy saved this way is
usually much less than the error due to discrete approximation,
so it is usually advisable not to try to get the absolute best
vertex positions.
<p>
There is one effect of hessian_normal that may be a little puzzling
at first. Many times a surface is known to have modes with zero
eigenvalue; translational or rotational modes, for example. Yet
no zero eigenvalues are reported. For example, with the cube
eigenvalues found above,
<pre> Eigencounts: 0 <, 0 ==, 193 >
1. 0.001687468263013
2. 0.001687468263013
3. 0.001687468263013
4. 0.2517282974725
5. 0.2517282974731
</pre>
one might expect 6 zero eigenvalues from three translational modes and
three rotational modes. But the rotational modes are eliminated
by the hessian_normal restriction. The three translational modes
have eigenvalue near 0, but not exactly 0, since normal motion can
approximate the translation of a cube, but not exactly.
The effective translation results
from vertices moving in on one hemisphere and out on the other.
This distorts the triangulation, raising the energy, hence the
positive eigenvalue. This effect decreases with refinement.
<p>
There are times when the normal direction is not the direction
one wants. If one knows the perturbation direction, one can use
the <a href="datafile.htm#hessian_special_normal_vector" class="keyword">
hessian_special_normal</a> feature to use that instead of the default
normal. Beware that
hessian_special_normal also applies to the normal calculated by
the <a href="elements.htm#vertex_normal" class="keyword">vertex_normal</a>
attribute and the normal
used by regular <a href="commands.htm#vertex_average">vertex averaging</a>.
<hr><a id="eigenvectorvisualization"></a><h2>
10. Visualizing eigenvectors.</h2>
Naturally, you want to see what various modes look like. At present.
Evolver can do this through a subsidiary menu invoked by the
command "<a href="commands.htm#hessian_menu" class="keyword">hessian_menu</a>".
This is a menu I use to test out various Hessianrelated features.
<table>
<tr>
<td>Option</td><td>Action</td></tr>
<tr>
<td>
1 </td><td>Initialize Hessian.</td></tr>
<tr>
<td style="verticalalign:texttop">
Z </td><td> Do ritz calculation of multiple eigenpairs. This
prompts for a shift value. Pick a value near the
eigenvalues you want, preferably below them so the
shifted Hessian is positive definite. This also
prompts for a number of eigenpairs to do. Eigenvalue
approximations are printed as they converge. You
can stop the iterations with a keyboard interrupt.</td></tr>
<tr>
<td style="verticalalign:texttop" >
X </td><td>Pick which eigenvector you want to use for moving.
You can enter
the eigenvector by its number in the list from the Z option.
As a special bonus useful when there are multiple eigenvectors for an
eigenvalue, you can enter the vector as a linear combination of
eigenvectors, e.g. "0.4 v1 + 1.3 v2  2.13 v3".</td></tr>
<tr>
<td style="verticalalign:texttop">
4 </td><td> Move. This prompts you for a scale factor. 1 is a
reasonable first guess. If it is too big or too
small, option 7 restores the original surface so
you can try option 4 again.</td></tr>
<tr>
<td style="verticalalign:texttop" >
7 </td><td>Restore original surface. Do this unless you want
your surface to stay moved when you exit hessian_menu.
Can repeat cycle by doing option X again.</td></tr>
<tr>
<td style="verticalalign:texttop" >
q </td><td>Quit back to main prompt. The surface is not
automatically restored to its original state. Be
sure to do option 7 before quitting if you don't
want to keep the changes!</td></tr>
</table>
<p>
<b>Square example.</b> Let us see how to display the lower eigenmodes
for the square, pictured here: <p>
<img src="squaree1.gif" alt="square eigenmode 1">
<img src="squaree2.gif" alt="square eigenmode 2">
<img src="squaree3.gif" alt="square eigenmode 3">
<img src="squaree4.gif" alt="square eigenmode 4">
<img src="squaree5.gif" alt="square eigenmode 5">
<img src="squaree6.gif" alt="square eigenmode 6">
<img src="squaree7.gif" alt="square eigenmode 7">
<img src="squaree8.gif" alt="square eigenmode 8">
<br>
Run square.fe, display, and refine three times. Run "linear_metric".
Run "hessian_menu". At the hessian_menu prompt, do option 1 (which
calculates the Hessian matrix),
then option z (ritz command) and respond to the prompt for number
of eigenvectors with 10, then option x and pick eigenvector 1.
Then pick option 4 and Step size 1. You should see the mode displayed
in the graphics window. Your image may look the opposite of the first
one pictured above, since the eigenvector comes out of a random initial
choice and thus may point in the opposite direction of what I got.
<p>After admiring the mode a while, do option 7 to restore the surface.
Then do option 4 again with Step size 1 to see the opposite mode.
Do option 7 to get back to the orignal. By doing the cycle of options
x, 4, and 7 you should be able to look at all the eigenmodes found
by the ritz command.
<p> The eigenvector components used to move vertices are stored
in the v_velocity vertex vector attribute, so you may access them
after exiting hessian_menu, for example <pre>
set vertex __x __x + 0.5*v_velocity
</pre>
<b>Raggedlooking eigenvectors:</b> If you have a volume constraint
and make a large move to show an eigenvector, you may get a bumpy
shape instead of the smooth shape you were expecting. This is because
the surface is projected to the volume constraint after moving in
the eigenvector direction. The eigenvector itself is smooth, but
the projection back to the target volume uses the volume gradient,
which depends very much on the triangulation and need not be smooth.
To see the eigenvector without volume projection, from a command
prompt you can do <pre>
set vertex __x __x + 0.5*v_velocity
</pre>
<p>
Hessian_menu can do sundry other things, mostly for my use in
debugging. Wouldbe gurus can
consult the <a href="commands.htm#hessian_menu" class="keyword">hessian_menu</a>
command documentation for a smorgasbord of more things to do with the Hessian.
<hr><a id="hessianiteration"></a>
<h2>11. Hessian iteration.</h2>
Or how to converge really really fast.
<p>
Suppose we assume the quadratic approximation to the surface energy is a
good one. Then we can find an equilibrium point by solving for
motion Y that makes the energy gradient zero. Recalling that G and
H depend only on X, the energy gradient as a function of Y is
<pre>
grad E = G<sup>T</sup> + Y<sup>T</sup> H.
</pre>
So we want G<sup>T</sup> + Y<sup>T</sup> H = 0, or transposing,
<pre>
G + H Y = 0.
</pre>
Solving for Y gives
<pre>
Y =  H<sup>1</sup> G.
</pre>
This is actually the NewtonRaphson Method applied to the gradient.
The Evolver's "<a href="commands.htm#hessiancommand" class="keyword">hessian</a>"
command does this calculation and
motion. It works best when the surface is near an equilibrium
point. It doesn't matter if the equilibrium is a minimum,
a saddle, or a maximum. However, nearness does matter. Remember
we are dealing with thousands of variables, and you don't have
to be very far away from an equilibrium for the equilibrium
to not be within the scope of the quadratic approximation.
When it does work, <code>hessian</code> will converge very quickly,
4 or 5 iterations at most.
<p>
Example: Start with the cube datafile cube.fe.
Evolve with "<code class="shortcode">r; g 5; r; g 5;</code>". This gets
very close to an equilibrium. Doing <code>hessian</code> a couple of times
gets to the equilibrium:
<pre>
Enter command: hessian
1. area: 4.85807791572284 energy: 4.85807791572284
Enter command: hessian
1. area: 4.85807791158432 energy: 4.85807791158432
Enter command: hessian
1. area: 4.85807791158431 energy: 4.85807791158431
</pre>
So Hessian iteration converged in two steps. Furthermore, this is
a local minimum rather than a saddle point, since Evolver did not
complain about a nonpositive definite Hessian.
<p>NOTE: The <code>hessian</code> command will work with indefinite Hessians,
as long as there are no eigenvalues too close to zero.
The warning about nonpositive definite is for your information;
it is does not mean <code>hessian</code> has failed.
<p>
<hr> <a id="hessian_seek"></a>
<h2> 12. Hessian_seek.</h2>
Even when the Hessian is positive definite, a hessian iteration
may blow up since the surface is not near enough to the minimum
for the Hessian approximation to be valid. For this circumstance,
the <code>hessian_seek</code> command does the same calculation as the <code>hessian</code>
command, except it does a line search along the direction of motion
for the minimum energy instead of jumping directly to the supposed
equilibrium. Basically, it does the same thing as the '<code>g</code>' command in
optimizing mode, except using the hessian solution instead of the
gradient.
<p>
<b>Cube example.</b> Run cube.fe, and do "<code class="shortcode">g; r; hessian_seek"</code>.
<pre>Enter command: g
0. area: 5.13907252918614 energy: 5.13907252918614 scale: 0.185675
Enter command: r
Vertices: 50 Edges: 144 Facets: 96 Bodies: 1 Facetedges: 288
Element memory: 40312
Total data memory: 299082 bytes.
Enter command: hessian_seek
1. area: 4.91067149153091 energy: 4.91067149153091 scale: 0.757800
</pre>
This is kind of a trivial example, since <code>hessian</code> doesn't blow up here.
But in general, when in doubt that hessian will work, use <code>hessian_seek</code>.
If <code>hessian_seek</code> reports a scale near 1, then it is probably safe to
use <code>hessian</code> to get the last few decimal places of convergence.
<code>Hessian_seek</code> cannot be as accurate as <code>hessian</code> at final convergence,
since <code>hessian_seek</code> uses differences of energies, which are very
small near the minimum.
<p>
<code>Hessian_seek</code> is safe to use even when the Hessian is not positive
definite. I have sometimes found it surprisingly effective even
when the surface is nowhere near a minimum.
<hr><a id="saddleexample"></a>
<h2> 13. Saddle.</h2>
Sometimes your surface may reach a saddle point of energy where gradient
descent becomes very very slow or even useless. Or you may be wanting
to use hessian, but there are pesky negative eigenvalues. The obvious
thing to do is pick a negative eigenvalue and move in the eigenvector
direction. You could do it by hand with <code>hessian_menu</code>, but the
<code>saddle</code>
command packages it all nicely. It finds the most negative eigenvector
and does a line search in that direction for the lowest energy.
To see an example, run catbody.fe with the following evolution:
<p> <img src="catbodysaddle.gif" alt="catenoid">
<img src="catbody.gif" alt="catenoid">
<pre> u
body[1].target := 4
g 5
r
g 5
V
V
r
g 10
eigenprobe 0
saddle
</pre>
Further evolution after <code>saddle</code> is necessary to find the lowest
energy surface; <code>saddle</code> just gives an initial push.
It is not necessary to run "<code class="shortcode">eigenprobe 0</code>" first;
if <code>saddle</code> finds no
negative eigenvalues, it says so and exits without moving the surface.
<p>
Saddle is safe to use even when nowhere near an equilibrium point.
<hr><a id="unstablesurfaces"></a><h2>
14. Detecting the onset of instability and evolving unstable surfaces.</h2>
The Hessian features can be used to detect the onset of instability
as some parameter changes, and even evolve unstable equilibrium surfaces.
<p>
Instability detection is done by watching eigenvalues with the
<a href="commands.htm#ritz" class="keyword">ritz</a>
command. As an example, consider a ring of liquid outside a cylinder, with
the volume increasing until the symmetric ring becomes unstable. This
is set up in the datafile catbody.fe, which is just cat.fe with a body
defined from the facets. Run catbody.fe with this initial evolution:
<pre> u
g 5
r
g 5
body[1].target := 2
g 5
r
body[1].target := 3
g 5
hessian
hessian
linear_metric
ritz(0,5)
</pre>This gives eigenvalues<pre>
Eigencounts: 0 <, 0 ==, 167 >
1. 0.398411128930840
2. 0.398411128930842
3. 1.905446082321839
4. 1.905446082321843
5. 4.4342055632012
</pre>
Note we are still in a stable, positive definite situation, but the
lowest eigenvalues are near enough to zero that we need to take care
in increasing the volume. Try an increment of 0.1:
<pre> body[1].target += 0.1
g 5
hessian
hessian
ritz(0,5)
Eigencounts: 0 <, 0 ==, 167 >
1. 0.287925880010193
2. 0.287925880010195
3. 1.775425717998147
4. 1.775425717998151
5. 4.2705109310529
</pre>
A little linear interpolation suggests try an increment of 0.3:
<pre> body[1].target += 0.3
g 5
hessian
hessian
hessian
ritz(0,5)
Eigencounts: 0 <, 0 ==, 167 >
1. 0.001364051154697
2. 0.0013640511547
3. 1.4344757227809
4. 1.4344757227809
5. 3.8350719808531
</pre>
So we are now very very close to the critical volume. In view of
the coarse triangulation here, it is probably not worth the trouble
to narrow down the critical volume further, but rather refine
and repeat the process. But for now keep this surface running
for the next paragraph.
<p>
<b>Evolving into unstable territory</b>
Typically these are surfaces with just a few unstable modes.
The idea is to get close to the desired equilibrium and use
<code>hessian</code> to reach it.
Regular '<code>g</code>' gradient descent
iteration should not be used. To change the surface,
i.e. to follow an equilibrium curve through a phase diagram,
make small changes to the control parameter and use a couple of
hessian iterations to reconverge each time. Particular care
is needed near bifurcation points, because of the several
equilibrium curves that meet. When approaching a bifurcation
point, try to jump over it so there is no ambiguity as to
which curve you are following. The approach to a bifurcation point
can be detected by watching eigenvalues. An eigenvalue crosses 0
when a surface introduces new modes of instability.
<p>
Example: Catbody.fe, continued. With catbody.fe in the nearlycritical
state found above, increase the body volume by steps of .1 and run
<code>hessian</code> a couple of times each step:
<pre> body[1].target += .1
hessian
hessian
hessian
hessian
body[1].target += .1
hessian
hessian
hessian
hessian
body[1].target += .1
hessian
hessian
hessian
hessian
</pre>
So <code>hessian</code> alone is enough to evolve with, as long as you stay
near enough to the equilibrium.
<p>
Other methods for evolving unstable surfaces:
<ul>
<li>Using symmetry. If the unstable surface is more symmetric than the
stable surfaces, then enforcing symmetry can remove the unstable modes.
For example, a surface of revolution could be retricted to just
a 90 degree wedge between two perpendicular mirror planes (levelset
constraints), with 90 degree contact angles on the planes.
<li>Using volume constraints. Recall that in general every constraint
removes one degree of freedom in the configuration space. Hence
a volume constraint has the potential to remove one unstable mode.
For example, unstable catenoids can be made stable by adding a
volume constraint and adjusting the volume until the pressure is 0.
</ul>
<hr><a id="eigenvalueaccuracy"></a>
<h2> 15. Eigenvalue Accuracy.</h2>
For the hardcore Evolver guru.
<p>
The question here is how well the eigenvalues of the discrete surface
approximate the eigenvalues of the smooth surface. This is significant
when deciding the stability of a surface. I present some comparisons
in cases where an analytic result is possible. All are done in
<a href="#hessiannormalmode" class="keyword">hessian_normal</a> mode with
<a href="#hessianmetric" class="keyword">linear_metric</a>.
<p>
Any general theorem on accuracy is difficult, since any of
these situations may arise:
<ul>
<li> The discrete surface exists, but the smooth surface does not.
<li> The smooth surface exists, but the discrete surface does not.
<li> The smooth and discrete surfaces exist, but their eigenvalues
are not close due to proximity to a bifurcation point.
<li> The discrete surface is a good approximation to the smooth
surface, and their eigenvalues are close. In general, linear model
eigenvalues will have error on the order h<sup>2</sup>, where h is a typical
edge length, and quadratic model eigenvalues will have error h<sup>4</sup>.
</ul>
<p>
<b>Example:</b> Flat square membrane with fixed sides of length pi. <br>
Smooth eigenvalues: n<sup>2</sup> + m<sup>2</sup>,
or 2,5,5,8,10,10,13,13,17,17,18,...
<p>
<pre>
64 facets:
linear_metric_mix:=0 linear_metric_mix:=.5 linear_metric_mix:=1
1. 1.977340485013809 2.036549240354332 2.098934776481970
2. 4.496631115530167 4.959720306550873 5.522143605551107
3. 4.496631115530168 4.959720306550873 5.522143605551109
4. 6.484555753109618 7.764634143334637 9.618809107463317
5. 8.383838160213188 10.098798069316681 12.451997407229225
6. 8.958685299442784 10.499717069102577 12.677342602918362
7. 9.459695221455723 12.274525789880890 17.295660791788656
8. 9.459695221455725 12.274525789880887 17.295660791788663
9. 11.5556960351898 15.7679635512509 23.585015056138165
10. 12.9691115062193 16.7214142904454 23.5850152363232
256 facets:
linear_metric_mix:=0 linear_metric_mix:=.5 linear_metric_mix:=1
1. 1.994813709598124 2.009974216370146 2.025331529636075
2. 4.869774462491779 4.999499042685448 5.135641457408189
3. 4.869774462491777 4.999499042685451 5.135641457408191
4. 7.597129628414264 7.985384943789075 8.411811839672165
5. 9.571559289947587 10.090156419079083 10.639318311067063
6. 9.759745949169531 10.186524915471141 10.665025703718413
7. 12.026114091326091 13.008227978344539 14.149533399632906
8. 12.026114091326104 13.008227978344539 14.149533399632912
9. 15.899097189772919 17.242998229925817 18.8011199268796
10. 15.8990975402264 17.2429983900303 18.8011200311777
1024 facets:
linear_metric_mix:=0 linear_metric_mix:=.5 linear_metric_mix:=1
1. 1.998755557957786 2.002568891165917 2.006394465437547
2. 4.967163232244397 5.0005039961225 5.0342462883517
3. 4.967163232244401 5.0005039961225 5.0342462883517
4. 7.897718646133286 7.9990889953386 8.1028624365882
5. 9.891301374223204 10.0268080202750 10.1637119963890
6. 9.941758861492074 10.0519390576227 10.1658353297015
7. 12.750608847206168 13.0141591049184 13.2878054981609
8. 12.750608847206173 13.0141591049184 13.2878054981609
9. 16.718575011911319 17.0839068189583 17.4625185925160
10. 16.7185751321348 17.0839069326949 17.4625186893608
</pre>
A curious fact here is that eigenvalues 5 and 6 are identical in
the smooth limit, but are split in the discrete approximation.
This is due to the fact that the 2dimension eigenspace of the
eigenvalue 10 can have a basis of two perturbations that each
have the symmetry of the square, but are differently affected by
the discretization.
<p>
<b>Example:</b> Sphere of radius 1 with fixed volume. <br>
Analytic eigenvalues: (n+2)(n1), n=1,2,3,...<br>
with multiplicities: 0,0,0,4,4,4,4,4,10,10,10,10,10,10,10
<pre>
Linear mode, linear_metric_mix:=.5.
24 facets 96 facets 384 facets 1536 facets
1. 0.3064952760319 0.1013854786956 0.0288140848394 0.0078006105505
2. 0.3064952760319 0.1013854786956 0.0288140848394 0.0078006105505
3. 0.3064952760319 0.1013854786956 0.0288140848394 0.0078006105505
4. 2.7132437122162 3.9199431944791 4.0054074145696 4.0036000699357
5. 2.7132437122162 3.9199431944791 4.0054074145696 4.0036000699357
6. 2.7132437122162 3.9199431944791 4.0054074145696 4.0036000699357
7. 3.6863290283837 4.3377418342267 4.1193008989205 4.0347069118257
8. 4.7902142880145 4.3377418342267 4.1193008989205 4.0347069118257
9. 4.7902142880145 9.0031399793203 9.9026247196089 9.9870390309740
10. 6.7380098215325 9.7223042306475 10.0981057119891 10.0387404054572
11. 6.7380098215325 9.7223042306545 10.0981057120367 10.0387404054607
12. 6.7380098215325 9.7223042307334 10.0981057121432 10.0387404055516
13. 8.6898276285107 9.9763109502799 10.1298354477703 10.0466252566958
In quadratic mode (net 4 times as many vertices per facet)
24 facets 96 facets 384 facets 1536 facets
1. 0.0311952242811 0.0025690857791 0.0002802874477 0.0000358409262
2. 0.0311952242812 0.0025690857791 0.0002802874477 0.0000358409262
3. 0.0311952242812 0.0025690857791 0.0002802874477 0.0000358409262
4. 4.2142037235384 4.0160946848738 4.0009978658738 4.0000515986034
5. 4.2142037235384 4.0160946848738 4.0009978658738 4.0000515986034
6. 4.2564592100813 4.0222815310390 4.0014892600742 4.0000883321792
7. 4.2564592100813 4.0222815310390 4.0014892600742 4.0000883321792
8. 4.2564592100813 4.0222815310390 4.0014892600742 4.0000883321792
9. 9.9900660130172 10.1007993043211 10.0076507048408 10.0004402861468
10. 9.9900660130172 10.1007993043271 10.0076507049338 10.0004402861545
11. 9.9900660130173 10.1007993047758 10.0076507050506 10.0004402861829
12. 12.1019587923328 10.1262981041797 10.0089922470136 10.0005516796159
13. 12.1019587923350 10.1262981042009 10.0089922470510 10.0005516796196
14. 12.1019587927495 10.1262981044110 10.0089922470904 10.0005516797052
15. 12.6934178610912 10.1478397915001 10.0106695599620 10.0007050467653
</pre>
<hr><a id="hessiancaveats"></a><h2>
16. Caveats and warnings.</h2>
When not near enough to an equilibrium, incautious use of
<a href="commands.htm#hessiancommand" class="keyword">hessian</a>
can wreck a surface. If you don't know what's going to happen,
save the surface first. Or use hessian_seek.
Or use <a href="commands.htm#hessian_menu" class="keyword">
hessian_menu</a>, where you can restore
the surface with option 7. Or set <a href="toggle.htm#check_increase" class="keyword">
check_increase</a>, which will
abort a hessian iteration that would increase energy. But remember
sometimes hessian should increase energy, for example to satisfy
a volume constraint or reach an unstable equilibrium.
<p>
Not all energies have builtin Hessians. If Evolver complains about
lacking a Hessian for a particular energy, you will either have
to forego <code>hessian</code> or find a way to phrase your problem in terms
of an energy with a Hessian.
<p>
There are three methods of sparse matrix factoring built into Evolver:
YSMP (Yale Sparse Matrix Package), my own minimal degree algorithm, and
<a href="install.htm#metisinstall">METIS</a> recursive bisection.
The <a href="toggle.htm#ysmp" class="keyword">ysmp</a> and
<a href="toggle.htm#metis_factor" class="keyword">
metis_factor</a> toggles can be used to control which is active.
<code>Metis</code> is the best overall, particularly on large surfaces, but requires
a separate download and compilation of the METIS library (but it's easy;
see the <a href="install.htm#metisinstall">installation</a> instructions).
<hr>
<a href="#hessiantutorial">back to Hessian top</a>
<hr>
<a href="evolver.htm#doctop">Back to top of Surface Evolver documentation.</a>
<a href="index.htm">Index.</a>
</body>
</html>
