1 | ### eigen =====================================================================
|
2 |
|
3 | Tags
|
4 | ----
|
5 | scripting, JS, internal, treenode, general concept
|
6 |
|
7 | Parameters
|
8 | ----------
|
9 | m
|
10 |
|
11 | General description
|
12 | -------------------
|
13 | Compute eigenvalues and eigenvectors. Matrix m must be both numerical and symmetric.
|
14 | The eigenval function returns a matrix with the eigenvalues along the diagonal.
|
15 | The eigenvec function returns a matrix with the eigenvectors arranged as row vectors.
|
16 | The eigen function does not return anything but stores the eigenvalue matrix in D
|
17 | and the eigenvector matrix in Q.
|
18 |
|
19 | Input: stack[tos - 1] symmetric matrix
|
20 |
|
21 | Output: D diagnonal matrix
|
22 | Q eigenvector matrix
|
23 |
|
24 | D and Q have the property that
|
25 |
|
26 | A == dot(transpose(Q),D,Q)
|
27 |
|
28 | where A is the original matrix.
|
29 |
|
30 | The eigenvalues are on the diagonal of D.
|
31 | The eigenvectors are row vectors in Q.
|
32 |
|
33 | The eigenvalue relation:
|
34 |
|
35 | A X = lambda X
|
36 |
|
37 | can be checked as follows:
|
38 |
|
39 | lambda = D[1,1]
|
40 | X = Q[1]
|
41 | dot(A,X) - lambda X
|
42 |
|
43 | Example 1. Check the relation AX = lambda X where lambda is an eigenvalue and X is the associated eigenvector.
|
44 |
|
45 | Enter:
|
46 |
|
47 | A = hilbert(3)
|
48 |
|
49 | eigen(A)
|
50 |
|
51 | lambda = D[1,1]
|
52 |
|
53 | X = Q[1]
|
54 |
|
55 | dot(A,X) - lambda X
|
56 |
|
57 | Result:
|
58 |
|
59 | -1.16435e-14
|
60 |
|
61 | -6.46705e-15
|
62 |
|
63 | -4.55191e-15
|
64 |
|
65 | Example 2: Check the relation A = QTDQ.
|
66 |
|
67 | Enter:
|
68 |
|
69 | A - dot(transpose(Q),D,Q)
|
70 |
|
71 | Result:
|
72 |
|
73 | 6.27365e-12 -1.58236e-11 1.81902e-11
|
74 |
|
75 | -1.58236e-11 -1.95365e-11 2.56514e-12
|
76 |
|
77 | 1.81902e-11 2.56514e-12 1.32627e-11
|
78 |
|
79 | ###
|
80 |
|
81 |
|
82 |
|
83 | #define D(i, j) yydd[EIG_N * (i) + (j)]
|
84 | #define Q(i, j) yyqq[EIG_N * (i) + (j)]
|
85 |
|
86 | EIG_N = 0
|
87 | EIG_yydd = []
|
88 | EIG_yyqq = []
|
89 |
|
90 | Eval_eigen = ->
|
91 | if (EIG_check_arg() == 0)
|
92 | stop("eigen: argument is not a square matrix")
|
93 |
|
94 | eigen(EIGEN)
|
95 |
|
96 | p1 = usr_symbol("D")
|
97 | set_binding(p1, p2)
|
98 |
|
99 | p1 = usr_symbol("Q")
|
100 | set_binding(p1, p3)
|
101 |
|
102 | push(symbol(NIL))
|
103 |
|
104 | ### eigenval =====================================================================
|
105 |
|
106 | Tags
|
107 | ----
|
108 | scripting, JS, internal, treenode, general concept
|
109 |
|
110 | Parameters
|
111 | ----------
|
112 | m
|
113 |
|
114 | General description
|
115 | -------------------
|
116 | Compute eigenvalues of m. See "eigen" for more info.
|
117 |
|
118 | ###
|
119 |
|
120 | Eval_eigenval = ->
|
121 | if (EIG_check_arg() == 0)
|
122 | push_symbol(EIGENVAL)
|
123 | push(p1)
|
124 | list(2)
|
125 | return
|
126 |
|
127 | eigen(EIGENVAL)
|
128 |
|
129 | push(p2)
|
130 |
|
131 | ### eigenvec =====================================================================
|
132 |
|
133 | Tags
|
134 | ----
|
135 | scripting, JS, internal, treenode, general concept
|
136 |
|
137 | Parameters
|
138 | ----------
|
139 | m
|
140 |
|
141 | General description
|
142 | -------------------
|
143 | Compute eigenvectors of m. See "eigen" for more info.
|
144 |
|
145 | ###
|
146 |
|
147 | Eval_eigenvec = ->
|
148 | if (EIG_check_arg() == 0)
|
149 | push_symbol(EIGENVEC)
|
150 | push(p1)
|
151 | list(2)
|
152 | return
|
153 |
|
154 | eigen(EIGENVEC)
|
155 |
|
156 | push(p3)
|
157 |
|
158 | EIG_check_arg = ->
|
159 | i = 0
|
160 | j = 0
|
161 |
|
162 | push(cadr(p1))
|
163 | Eval()
|
164 | yyfloat()
|
165 | Eval()
|
166 | p1 = pop()
|
167 |
|
168 | if (!istensor(p1))
|
169 | return 0
|
170 |
|
171 | if (p1.tensor.ndim != 2 || p1.tensor.dim[0] != p1.tensor.dim[1])
|
172 | stop("eigen: argument is not a square matrix")
|
173 |
|
174 | EIG_N = p1.tensor.dim[0]
|
175 |
|
176 | for i in [0...EIG_N]
|
177 | for j in [0...EIG_N]
|
178 | if (!isdouble(p1.tensor.elem[EIG_N * i + j]))
|
179 | stop("eigen: matrix is not numerical")
|
180 |
|
181 | for i in [0...(EIG_N-1)]
|
182 | for j in [(i+1)...EIG_N]
|
183 | if (Math.abs(p1.tensor.elem[EIG_N * i + j].d - p1.tensor.elem[EIG_N * j + i].d) > 1e-10)
|
184 | stop("eigen: matrix is not symmetrical")
|
185 |
|
186 | return 1
|
187 |
|
188 | #-----------------------------------------------------------------------------
|
189 | #
|
190 | # Input: p1 matrix
|
191 | #
|
192 | # Output: p2 eigenvalues
|
193 | #
|
194 | # p3 eigenvectors
|
195 | #
|
196 | #-----------------------------------------------------------------------------
|
197 |
|
198 | eigen = (op) ->
|
199 | i = 0
|
200 | j = 0
|
201 |
|
202 | # malloc working vars
|
203 |
|
204 | #EIG_yydd = (double *) malloc(n * n * sizeof (double))
|
205 | for i in [0...(EIG_N*EIG_N)]
|
206 | EIG_yydd[i] = 0.0
|
207 |
|
208 | #if (EIG_yydd == NULL)
|
209 | # stop("malloc failure")
|
210 |
|
211 | #EIG_yyqq = (double *) malloc(n * n * sizeof (double))
|
212 | for i in [0...(EIG_N*EIG_N)]
|
213 | EIG_yyqq[i] = 0.0
|
214 |
|
215 | #if (EIG_yyqq == NULL)
|
216 | # stop("malloc failure")
|
217 |
|
218 | # initialize D
|
219 |
|
220 | for i in [0...EIG_N]
|
221 | EIG_yydd[EIG_N * (i) + (i)] = p1.tensor.elem[EIG_N * i + i].d
|
222 | for j in [(i+1)...EIG_N]
|
223 | EIG_yydd[EIG_N * (i) + (j)] = p1.tensor.elem[EIG_N * i + j].d
|
224 | EIG_yydd[EIG_N * (j) + (i)] = p1.tensor.elem[EIG_N * i + j].d
|
225 |
|
226 | # initialize Q
|
227 |
|
228 | for i in [0...EIG_N]
|
229 | EIG_yyqq[EIG_N * (i) + (i)] = 1.0
|
230 | for j in [(i+1)...EIG_N]
|
231 | EIG_yyqq[EIG_N * (i) + (j)] = 0.0
|
232 | EIG_yyqq[EIG_N * (j) + (i)] = 0.0
|
233 |
|
234 | # step up to 100 times
|
235 |
|
236 | for i in [0...100]
|
237 | if (step() == 0)
|
238 | break
|
239 |
|
240 | if (i == 100)
|
241 | printstr("\nnote: eigen did not converge\n")
|
242 |
|
243 | # p2 = D
|
244 |
|
245 | if (op == EIGEN || op == EIGENVAL)
|
246 |
|
247 | push(p1)
|
248 | copy_tensor()
|
249 | p2 = pop()
|
250 |
|
251 | for i in [0...EIG_N]
|
252 | for j in [0...EIG_N]
|
253 | push_double(EIG_yydd[EIG_N * (i) + (j)])
|
254 | p2.tensor.elem[EIG_N * i + j] = pop()
|
255 |
|
256 | # p3 = Q
|
257 |
|
258 | if (op == EIGEN || op == EIGENVEC)
|
259 |
|
260 | push(p1)
|
261 | copy_tensor()
|
262 | p3 = pop()
|
263 |
|
264 | for i in [0...EIG_N]
|
265 | for j in [0...EIG_N]
|
266 | push_double(EIG_yyqq[EIG_N * (i) + (j)])
|
267 | p3.tensor.elem[EIG_N * i + j] = pop()
|
268 |
|
269 | # free working vars
|
270 |
|
271 |
|
272 | #-----------------------------------------------------------------------------
|
273 | #
|
274 | # Example: p = 1, q = 3
|
275 | #
|
276 | # c 0 s 0
|
277 | #
|
278 | # 0 1 0 0
|
279 | # G =
|
280 | # -s 0 c 0
|
281 | #
|
282 | # 0 0 0 1
|
283 | #
|
284 | # The effect of multiplying G times A is...
|
285 | #
|
286 | # row 1 of A = c (row 1 of A ) + s (row 3 of A )
|
287 | # n+1 n n
|
288 | #
|
289 | # row 3 of A = c (row 3 of A ) - s (row 1 of A )
|
290 | # n+1 n n
|
291 | #
|
292 | # In terms of components the overall effect is...
|
293 | #
|
294 | # row 1 = c row 1 + s row 3
|
295 | #
|
296 | # A[1,1] = c A[1,1] + s A[3,1]
|
297 | #
|
298 | # A[1,2] = c A[1,2] + s A[3,2]
|
299 | #
|
300 | # A[1,3] = c A[1,3] + s A[3,3]
|
301 | #
|
302 | # A[1,4] = c A[1,4] + s A[3,4]
|
303 | #
|
304 | # row 3 = c row 3 - s row 1
|
305 | #
|
306 | # A[3,1] = c A[3,1] - s A[1,1]
|
307 | #
|
308 | # A[3,2] = c A[3,2] - s A[1,2]
|
309 | #
|
310 | # A[3,3] = c A[3,3] - s A[1,3]
|
311 | #
|
312 | # A[3,4] = c A[3,4] - s A[1,4]
|
313 | #
|
314 | # T
|
315 | # The effect of multiplying A times G is...
|
316 | #
|
317 | # col 1 of A = c (col 1 of A ) + s (col 3 of A )
|
318 | # n+1 n n
|
319 | #
|
320 | # col 3 of A = c (col 3 of A ) - s (col 1 of A )
|
321 | # n+1 n n
|
322 | #
|
323 | # In terms of components the overall effect is...
|
324 | #
|
325 | # col 1 = c col 1 + s col 3
|
326 | #
|
327 | # A[1,1] = c A[1,1] + s A[1,3]
|
328 | #
|
329 | # A[2,1] = c A[2,1] + s A[2,3]
|
330 | #
|
331 | # A[3,1] = c A[3,1] + s A[3,3]
|
332 | #
|
333 | # A[4,1] = c A[4,1] + s A[4,3]
|
334 | #
|
335 | # col 3 = c col 3 - s col 1
|
336 | #
|
337 | # A[1,3] = c A[1,3] - s A[1,1]
|
338 | #
|
339 | # A[2,3] = c A[2,3] - s A[2,1]
|
340 | #
|
341 | # A[3,3] = c A[3,3] - s A[3,1]
|
342 | #
|
343 | # A[4,3] = c A[4,3] - s A[4,1]
|
344 | #
|
345 | # What we want to do is just compute the upper triangle of A since we
|
346 | # know the lower triangle is identical.
|
347 | #
|
348 | # In other words, we just want to update components A[i,j] where i < j.
|
349 | #
|
350 | #-----------------------------------------------------------------------------
|
351 | #
|
352 | # Example: p = 2, q = 5
|
353 | #
|
354 | # p q
|
355 | #
|
356 | # j=1 j=2 j=3 j=4 j=5 j=6
|
357 | #
|
358 | # i=1 . A[1,2] . . A[1,5] .
|
359 | #
|
360 | # p i=2 A[2,1] A[2,2] A[2,3] A[2,4] A[2,5] A[2,6]
|
361 | #
|
362 | # i=3 . A[3,2] . . A[3,5] .
|
363 | #
|
364 | # i=4 . A[4,2] . . A[4,5] .
|
365 | #
|
366 | # q i=5 A[5,1] A[5,2] A[5,3] A[5,4] A[5,5] A[5,6]
|
367 | #
|
368 | # i=6 . A[6,2] . . A[6,5] .
|
369 | #
|
370 | #-----------------------------------------------------------------------------
|
371 | #
|
372 | # This is what B = GA does:
|
373 | #
|
374 | # row 2 = c row 2 + s row 5
|
375 | #
|
376 | # B[2,1] = c * A[2,1] + s * A[5,1]
|
377 | # B[2,2] = c * A[2,2] + s * A[5,2]
|
378 | # B[2,3] = c * A[2,3] + s * A[5,3]
|
379 | # B[2,4] = c * A[2,4] + s * A[5,4]
|
380 | # B[2,5] = c * A[2,5] + s * A[5,5]
|
381 | # B[2,6] = c * A[2,6] + s * A[5,6]
|
382 | #
|
383 | # row 5 = c row 5 - s row 2
|
384 | #
|
385 | # B[5,1] = c * A[5,1] + s * A[2,1]
|
386 | # B[5,2] = c * A[5,2] + s * A[2,2]
|
387 | # B[5,3] = c * A[5,3] + s * A[2,3]
|
388 | # B[5,4] = c * A[5,4] + s * A[2,4]
|
389 | # B[5,5] = c * A[5,5] + s * A[2,5]
|
390 | # B[5,6] = c * A[5,6] + s * A[2,6]
|
391 | #
|
392 | # T
|
393 | # This is what BG does:
|
394 | #
|
395 | # col 2 = c col 2 + s col 5
|
396 | #
|
397 | # B[1,2] = c * A[1,2] + s * A[1,5]
|
398 | # B[2,2] = c * A[2,2] + s * A[2,5]
|
399 | # B[3,2] = c * A[3,2] + s * A[3,5]
|
400 | # B[4,2] = c * A[4,2] + s * A[4,5]
|
401 | # B[5,2] = c * A[5,2] + s * A[5,5]
|
402 | # B[6,2] = c * A[6,2] + s * A[6,5]
|
403 | #
|
404 | # col 5 = c col 5 - s col 2
|
405 | #
|
406 | # B[1,5] = c * A[1,5] - s * A[1,2]
|
407 | # B[2,5] = c * A[2,5] - s * A[2,2]
|
408 | # B[3,5] = c * A[3,5] - s * A[3,2]
|
409 | # B[4,5] = c * A[4,5] - s * A[4,2]
|
410 | # B[5,5] = c * A[5,5] - s * A[5,2]
|
411 | # B[6,5] = c * A[6,5] - s * A[6,2]
|
412 | #
|
413 | #-----------------------------------------------------------------------------
|
414 | #
|
415 | # Step 1: Just do upper triangle (i < j), B[2,5] = 0
|
416 | #
|
417 | # B[1,2] = c * A[1,2] + s * A[1,5]
|
418 | #
|
419 | # B[2,3] = c * A[2,3] + s * A[5,3]
|
420 | # B[2,4] = c * A[2,4] + s * A[5,4]
|
421 | # B[2,6] = c * A[2,6] + s * A[5,6]
|
422 | #
|
423 | # B[1,5] = c * A[1,5] - s * A[1,2]
|
424 | # B[3,5] = c * A[3,5] - s * A[3,2]
|
425 | # B[4,5] = c * A[4,5] - s * A[4,2]
|
426 | #
|
427 | # B[5,6] = c * A[5,6] + s * A[2,6]
|
428 | #
|
429 | #-----------------------------------------------------------------------------
|
430 | #
|
431 | # Step 2: Transpose where i > j since A[i,j] == A[j,i]
|
432 | #
|
433 | # B[1,2] = c * A[1,2] + s * A[1,5]
|
434 | #
|
435 | # B[2,3] = c * A[2,3] + s * A[3,5]
|
436 | # B[2,4] = c * A[2,4] + s * A[4,5]
|
437 | # B[2,6] = c * A[2,6] + s * A[5,6]
|
438 | #
|
439 | # B[1,5] = c * A[1,5] - s * A[1,2]
|
440 | # B[3,5] = c * A[3,5] - s * A[2,3]
|
441 | # B[4,5] = c * A[4,5] - s * A[2,4]
|
442 | #
|
443 | # B[5,6] = c * A[5,6] + s * A[2,6]
|
444 | #
|
445 | #-----------------------------------------------------------------------------
|
446 | #
|
447 | # Step 3: Same as above except reorder
|
448 | #
|
449 | # k < p (k = 1)
|
450 | #
|
451 | # A[1,2] = c * A[1,2] + s * A[1,5]
|
452 | # A[1,5] = c * A[1,5] - s * A[1,2]
|
453 | #
|
454 | # p < k < q (k = 3..4)
|
455 | #
|
456 | # A[2,3] = c * A[2,3] + s * A[3,5]
|
457 | # A[3,5] = c * A[3,5] - s * A[2,3]
|
458 | #
|
459 | # A[2,4] = c * A[2,4] + s * A[4,5]
|
460 | # A[4,5] = c * A[4,5] - s * A[2,4]
|
461 | #
|
462 | # q < k (k = 6)
|
463 | #
|
464 | # A[2,6] = c * A[2,6] + s * A[5,6]
|
465 | # A[5,6] = c * A[5,6] - s * A[2,6]
|
466 | #
|
467 | #-----------------------------------------------------------------------------
|
468 |
|
469 | step = ->
|
470 | i = 0
|
471 | j = 0
|
472 | count = 0
|
473 |
|
474 | # for each upper triangle "off-diagonal" component do step2
|
475 |
|
476 | for i in [0...(EIG_N-1)]
|
477 | for j in [(i+1)...EIG_N]
|
478 | if (EIG_yydd[EIG_N * (i) + (j)] != 0.0)
|
479 | step2(i, j)
|
480 | count++
|
481 |
|
482 | return count
|
483 |
|
484 | step2 = (p,q) ->
|
485 | k = 0
|
486 | t = 0.0
|
487 | theta = 0.0
|
488 | c = 0.0
|
489 | cc = 0.0
|
490 | s = 0.0
|
491 | ss = 0.0
|
492 |
|
493 | # compute c and s
|
494 |
|
495 | # from Numerical Recipes (except they have a_qq - a_pp)
|
496 |
|
497 | theta = 0.5 * (EIG_yydd[EIG_N * (p) + (p)] - EIG_yydd[EIG_N * (q) + (q)]) / EIG_yydd[EIG_N * (p) + (q)]
|
498 |
|
499 | t = 1.0 / (Math.abs(theta) + Math.sqrt(theta * theta + 1.0))
|
500 |
|
501 | if (theta < 0.0)
|
502 | t = -t
|
503 |
|
504 | c = 1.0 / Math.sqrt(t * t + 1.0)
|
505 |
|
506 | s = t * c
|
507 |
|
508 | # D = GD
|
509 |
|
510 | # which means "add rows"
|
511 |
|
512 | for k in [0...EIG_N]
|
513 | cc = EIG_yydd[EIG_N * (p) + (k)]
|
514 | ss = EIG_yydd[EIG_N * (q) + (k)]
|
515 | EIG_yydd[EIG_N * (p) + (k)] = c * cc + s * ss
|
516 | EIG_yydd[EIG_N * (q) + (k)] = c * ss - s * cc
|
517 |
|
518 | # D = D transpose(G)
|
519 |
|
520 | # which means "add columns"
|
521 |
|
522 | for k in [0...EIG_N]
|
523 | cc = EIG_yydd[EIG_N * (k) + (p)]
|
524 | ss = EIG_yydd[EIG_N * (k) + (q)]
|
525 | EIG_yydd[EIG_N * (k) + (p)] = c * cc + s * ss
|
526 | EIG_yydd[EIG_N * (k) + (q)] = c * ss - s * cc
|
527 |
|
528 | # Q = GQ
|
529 |
|
530 | # which means "add rows"
|
531 |
|
532 | for k in [0...EIG_N]
|
533 | cc = EIG_yyqq[EIG_N * (p) + (k)]
|
534 | ss = EIG_yyqq[EIG_N * (q) + (k)]
|
535 | EIG_yyqq[EIG_N * (p) + (k)] = c * cc + s * ss
|
536 | EIG_yyqq[EIG_N * (q) + (k)] = c * ss - s * cc
|
537 |
|
538 | EIG_yydd[EIG_N * (p) + (q)] = 0.0
|
539 | EIG_yydd[EIG_N * (q) + (p)] = 0.0
|
540 |
|
541 |
|