See edit.
Sorry for the noise then. I forgot it was the preview that you were showing.
1 Like
About the eig() function and eigen command:
- This works only for symmetric matrices (real eigenvalues).
- It uses SVD decomposition to do the computation
- The SVD routine has been mainly converted from the version proposed in the Numerical Recipes (in C).
- This routine is not parallelized (it seems to be very hard to parallelize efficiently BTW).
I must say Iâve never been that disapointed by the performance of it. Computing all eigenvalues of a large (general) matrix is a computationally-expensive task. Apart from very specialized math toolboxes (e.g. Octave or Matlab) that probably does a better job than GâMIC, I really doubt youâll get that much speed improvements with other generic image processing toolboxes (Iâm not even sure how many IP toolboxes propose this feature 
).
EDIT: It takes approx. 30 seconds here for a 1024x1024 matrix, and 3 seconds for a 512x512 one. If I remember well, the overall algorithm complexity is in O(N^3) where the matrix dimension is N\times N.
1 Like
Sooo almost a year of computation time for the well with matrix N = 300
1 Like
By the way, looking at the different matrices posted there : they are non general matrices, but band matrices, and it probably exists a faster algorithm to compute eigenvalues of such particular matrices than considering the general case.
(An improved divide-and-conquer algorithm for the banded matrices with narrow bandwidths - ScienceDirect maybe?).
In any case, GâMIC doesnât have a native command for this (but this is probably possible to write a custom command to handle this case).
1 Like
Is that why eigen didnât work out here? If it is, I think I will hold off this project then.
I donât know what you mean by âdidnât work outâ. Are you taking about the long processing time, or do you get inaccurate results ? I didnât follow the thread in details.
What I can say is that toolboxes for matrix computation may first analyze the type of matrix to process, and then, choose the âbestâ algorithm specifically tailored for a specific class of matrices.
Then, the user may feel computing the result is simple enough, but usually, itâs not. Itâs just that the math toolbox does the complex work for you 
Basically, I didnât get the result similar to the python result. I do think that I did followed the python script well going off the description. It is possible that I donât know how to use eigen as this is a first for me.
1 Like
Another bus trip to Schrodingerâs Well and - hey! - this time I got great pictures! Helps to bring the right camera, lens, filter, tripodâŚ
Time to put some concerns to bed:
-
Instability: Due to non-nonsensical seed map. A seedmap is a binary-value cross-sectional view of the well. Make a circle, square, warp it any way you like. I think it could even be genus torus-n for n, a bunch of holes. Once I had a decent seed map, eigs didnât crash anymore. Even at N=300. Processing time at N=300 was eight minutes in the guts of
eigs()- not optimized, so computation was in one core at a modest 2.2 GHz. Your mileage may vary. - Wrong Hamiltonian: Simple bug. Iâd modified @weightt_an code somewhat and dropped a wrench in the middle of it. @weightt_an picture helped straighten that out. Always keep tools not in use latched to the beltâŚ
- Output: There is a bit of post-processing that goes on where @weightt_an sorts the resulting eigenvalues/vectors for making an optimal choice. Looks straightforward to translate. Thoughts on that in the next post.
In raw terms, the output data are single-channel, ranging from ~0.01 to ~0.3. My post-processing consisted of renormalizing to 0,255, a slight gamma correction to favor the dark end - shifting that nearer to middle-gray, and translating via -map (used ârainbowâ).
User interface: @Reptorianâs department, but methinks it would be roughly (a) A binary image that gives a well cross-section. (b) A percentage hint for fast/low-resolution or slow/high-resolution. My snapshots are N=300, which internally calls for a 90000,90000,1,1 to hold the Hamiltonian image. The other significant choice is "How many energy levels (eigen values) you want to see? Methinks that ought to be internally fixed to 90-130, at it would be a subtle thing to convey to a casual user in any meaningful way.
N can be derived from the quality percentage the user gives, with 1% corresponding in some such way to N=30. The input image should be padded out to square and scaled so that N=300 means a re-scaled, padded input image of 300,300,1,1
Havenât hooked in GâMICâs eig() as a substitute for scipyâs eigs() but am now fairly confident that it will work for some useful value of N, maybe even N=300. Thank you @David_Tschumperle for the bit on eigâs implementation: know the Recipes in C version well, having stolen it for various purposes here and again. Probably next will be for me to walk through @Reptorian implementation, which looks generally right by eyeball inspection. Stay tunedâŚ
1 Like
FYI, eig() and eigen sorts the eigenvalues in decreasing order.
1 Like
Yes. I donât think âchoosing optimalâ is of great concern insofar as selecting from eigenvalues/vectors. I think the only thing between now and a wrap is giving GâMICâs -eigen a test drive, mainly to determine at what value N processing becomes unbearable. We may be pleasantly surprised there. You can fold Hamiltonians along the diagonal and the matching coefficients off-diagonal are equal, so we are not feeding -eigen anything strange.
One UI thing I neglected to mention is that there is a 1 â many relationship between input well cross-sections and output energy-level images. My 20 snapshots above all derived from one well cross-section and the eigenvalue/vector pairs coming out of one dive into eigs():
So @Reptorian can have a bunch of gray-scales to shift around - I look forward to what might be done with that.
1 Like
Another bus trip to Schrodingerâs Well! AndâŚ
âŚI fell in.
TL;DR: You donât want to do this. Seriously. Get a dog. Take walks in the park. Enjoy life.
Long Story: Generated the Hamiltonian matrix using @weightt_an script, but this time saved the H. array as a GâMIC .cimg for the purpose of loading it into a pipeline and running -eigen on the results.
Good news: -eigen solves for eigen values and vectors and the probability distributions do not differ from those generated by scipy.sparse.linalg.eigs() in @weightt_anâs script.
Bad news Lord⌠but -eigen does take its own sweet time doing it. For N=300, Iâll be finished writing tutorials for all 952+ GâMIC commands by the time it is done.
Timed by -tic -toc bracketing -eigen. Used a circular well seed image. N sizes, H. matrix image size, run times for -eigen in seconds:
N=30 - 900,900,1,1 - 7.1
N=40 - 1600,1600,1,1 - 107.5
N=50 - 2500,2500,1,1 - 736.
N=60 - 3600,3600,1,1 - 2710.
I think ten minutes is probably past the pain threshold of most people, (N=50) and the resulting images of probability distribution are tiny: 60x60 pixels on the last run, which was close to an hour. I decided not to sort or select the resulting vectors, but just squared each vector up into an image, normalized to 0,255 and packed the lot - 3600 vectors in the last run - and rolled them into a video clip. Itâs interesting, at an academic level, to see all the vectors and the patterns they make when squared up and normalized, paging through them in solution sequence. Have a look:
But the experience is like watching television at the 1939 New York Worldâs Fair. (before my time, but met people who were there.) The video, for me, was a proof of concept that I was getting images similar to what @weightt_an script generated, so I wasnât in left field, computationally.
The difficulty with -eigen compared to its scipy cousin centers primarily on being able to set the number of returned eigenvalue/eigenvector pairs, with the scipy implementation, which seems to short-circuit the computation time. The more eigen value/vector pairs you ask for with the scipy implementation, the longer running time goes, and scipy documentation - flat out - tells you that the implementation will not deliver all eigenvalue/vector pairs before the heat death of the Universe. The python script generally asked for about a 100 pairs, a small subset of the 900,000 that could be generated from an N=300 Hamiltonian.
In contrast, -eigen gives all eigenvalue/vector pairs, and youâll wait for all of them, even if you donât want all of them.
For interest, this is how I dumped from Python to a .cimg
Hamiltonian = diags(diagsV, diagsK, format = 'dia')
print('Hamiltonian values done')
if dumpham:
fn = 'Ham_{0:05d}.cimg'.format(No_points)
print('Dumping the Hamiltonian to {0:s}'.format(fn))
rwcimg.writecimg(fn, Hamiltonian.toarray())
e_values, e_vec = eigs(Hamiltonian, k = Elevels, ncv=4*Elevels, tol = mytol )
rwcimg is a module I wrote some time ago to write numpy arrays to disk as .cimg files. It is not much more than a copy routine. There is nothing strange or complicated about .cimg format. That -eigen reproduced the same probability distributions as @weightt_an script - just small - leads me to think that no artifacts were introduced there.
This is how I harnessed eigen:
gmic \
Ham_03600.cimg \
-tic \
-eigen. \
-toc \
-rm[0] \
-rotate. 90 \
-repeat {w} \
-shared. {$>},{$>},0,0 \
{sqrt(w)},{sqrt(w)},1,1 \
-fill. [-2] \
-remove[-2] \
sqr. \
n. 0,255 \
move. {$>} \
-done \
-remove. \
-output wellstack.mp4,16
-eigen solves the Hamiltonian, delivering two images: a 1,N*N,1,1 column vector of eigen values. These were removed from the image list; the probability distribution is actually imaged from the eigen vectors alone. These vectors come in the form of a N*N,N*N,1,1 image where each vector is an image column. I rotated the image 90 degrees so that the column vectors became row vectors - and occupying contiguous memory. That way, I could share out the rows one at a time and use the row to fill a square N,N,1,1 image, squaring and normalizing the vector components to a 0,255 range. Exiting the repeat loop, I dropped the eigen vector image and wrote out the N*N squared up images to a video.
So. The probability distributions are visually interesting, even compelling, but the road to get to them is a long one. The computational approach is likely naĂŻve and clever people will have to come to the fore with optimizations. Easier said then done. Ah, @Reptorian. Methinks you may be taking up another project. Sorry.
3 Likes
There is something clear there : you are definitely not interested in all the eigenvalues of the large matrices you handle. That makes a big difference, because computing the N first eigenvalues can be done with specific (faster) algorithms.
In this context, the article: https://arxiv.org/pdf/0909.4061.pdf seems to be an interesting read. Iâll check to see if that can be adapted easily as a GâMIC function.
Yes, I am in the process of adding foreground extraction selection to Krita with guided filtering technique.
Also, eigen doesnât seem to give me the result I want with my gmic translation. I would like to know why.
From what I see, the computation of the Nth first eigenvalues could be speed up by a significant factor, by using the Lanczos Algorithm, which seems completely doable as a GâMIC custom command.
I donât have much time right now to try, but it could be an interesting addition to GâMIC, so letâs put this on my todo list.
2 Likes
In fact, back in the day when I was using a crappy seed image and eigs() was tossing exceptions at me, thatâs what the error message referenced:
ArpackError: ARPACK error 3: No shifts could be applied during a cycle of the Implicitly restarted Arnoldi iteration. One possibility is to increase the size of NCV relative to NEV.
So. Hint, hint. Also, everything scipy is open source, so in principle it should be straight foward to study eigs() implementation. (NCV and NEV are tuning parameters. Thereâs a link to eigs() documentation upstream from here)
No surprised then, they use Arnoldi iterations as well
Iâm experimenting right now with the simplest âpower methodâ which already seems to be quite fast, and which could be enough for our application here.
1 Like
Looks about right (by eyeball). But never actually tried your GâMIC version of diags() - first wanted to know what kind of N could be gotten out of -eigen without becoming suicidal. So I was working backwards, toward you, but never got to your script. One more bus ride to Schrodingerâs WellâŚ