Timeline for Approximating sum of entries of $\exp(A-B)$ for diagonal $A$ and rank-$1$ $B$?
Current License: CC BY-SA 4.0
41 events
| when toggle format | what | by | license | comment | |
|---|---|---|---|---|---|
| Apr 21, 2023 at 5:56 | comment | added | Yaroslav Bulatov | btw, I gave background on the underlying problem and referenced your implementation here | |
| Apr 20, 2023 at 5:23 | comment | added | Yaroslav Bulatov | BTW, I've been wrapping my head around this method, and it appears this method comes out as special case of the residue theorem. Depending on the contour you choose, you will get either this method, or the standard method based on eigenvalue decomposition | |
| Mar 30, 2023 at 15:59 | comment | added | fedja | @YaroslavBulatov " if you have ideas how to do the larger scale version, I'll be happy to post this as a separate question" Sure, but since you are looking for something better than $O(d)$, you'll have to explain clearly what assumptions about the data you are making since in this case one cannot even read the data entirely | |
| Mar 30, 2023 at 2:44 | comment | added | Yaroslav Bulatov | So the solution you provided above solves the problem in notebook accurately in $O(d)$ time which works great medium problems. For the largest problems $d$ can be trillions, so I'd need an approximate solution it in $O(1)$ time, if you have ideas how to do the larger scale version, I'll be happy to post this as a separate question | |
| Mar 30, 2023 at 2:33 | comment | added | Yaroslav Bulatov | Thanks for providing the code + approximation pointers, it was quite useful for validating hypotheses and I learned a lot. I actually don't need any higher accuracy -- the reason is that matrices $A$,$B$ are not known and have to be estimated from data, sometimes all you know is that eigenvalues follow power-law decay with approximately known constant $p$. The underlying problem behind this question is summarized here. I can solve the "diag" case with explicit formula, but not "diag+rank1" case | |
| Mar 24, 2023 at 21:53 | comment | added | fedja | @YaroslavBulatov I added a few words to the answer to explain why the computation of $\lambda$ is absolutely indispensable if you want the code to work correctly, what exactly you should verify so that you would know what happens for large $t$ and $d$, not just assume it, and why your assumptions about the error in your code without $\lambda$ were completely unjustified. As usual, feel free to ask questions, if anything is unclear :-) | |
| Mar 24, 2023 at 21:48 | history | edited | fedja | CC BY-SA 4.0 |
added 2023 characters in body
|
| Mar 24, 2023 at 15:11 | comment | added | fedja | @YaroslavBulatov Actually, you'll be way off even with $d=100, t=10^{16}$, which you can compute directly. Just run your code as it is with that $d$ increasing $t$ along powers of 2 (squaring a 100 by 100 matrix is a fast enough operation). You'll see something interesting :-) | |
| Mar 24, 2023 at 14:21 | comment | added | fedja | "for d=10^8 I don't have a way to see what the error is, so I'll assume it's small." @YaroslavBulatov You do have a way, if you know where to look, and that is where $\lambda$ and its precision comes into play. We are not in church: there is no need to accept anything on faith. I'll explain it later today. But first, you should incorporate the computation of $\lambda$ and reduction by $\lambda I$ into your code. Without it, you may be fine for $d,t=100$, but will be way off for $d=10^8, t=10^{16}$. | |
| Mar 24, 2023 at 5:09 | vote | accept | Yaroslav Bulatov | ||
| Mar 24, 2023 at 5:09 | comment | added | Yaroslav Bulatov | No, I was evaluating with d=100 because it takes O(d^3) to compare relative error against ground truth (explicit computation of matrix exponential). For d=10^8 I don't have a way to see what the error is, so I'll assume it's small. I see that your method is very similar to Section 3 of Gallopoulos with the difference that they use Chebychev approximation. They give coefficients in Table 5 but seems like overkill, this solution appears already to have high enough accuracy, thanks! | |
| Mar 24, 2023 at 1:49 | comment | added | fedja | @YaroslavBulatov Here are the guaranteed relative errors that you should have if you implement the code correctly for n=16 (8 roots) in the format $(d,E)$: $(10^1, 2\cdot 10^{-8}),(10^2, 2.2\cdot 10^{-7}), (10^3, 2.2*10^{-6} ), (10^4, 2.1\cdot 10^{-5} ), (10^5, 1.6*10^{-4}), (10^6, 10^{-4} ), (10^7, 0.005), (10^8, 0.014 ), (10^9, 0.038)$. That should be true for all $t>0$ for which $t$ times the error you've made in the approximation of $\lambda$ is negligible. (so if you want it to be correct for $t=10^{18}$, make sure that you $\lambda$ is correct up to $10^{-25}$, say. | |
| Mar 24, 2023 at 0:50 | comment | added | fedja | @YaroslavBulatov Also, notice that I'm not trying to approximate $e^{-tV}$ but $e^{-tV/2}$. Yes, it is possible to write the answer as $(P^{-1}(V)w,w)$ and then you'll just need to sum the entries instead of their squares, but that would require a higher degree of the polynomial to get the decent precision and I was trying to keep it as low as feasible. Want to discuss it all on Skype/Zoom in person? | |
| Mar 24, 2023 at 0:41 | comment | added | fedja | @YaroslavBulatov Erm... What are you testing? $d=2$ with $t\in[0,100]$? That is not what you need 8 roots for! Also, I do not see the preliminary reduction by $\lambda I$, which is essential if you want small relative error in the whole range of t for large d. You should test $d$ up to $10^8$ and $t$ running up to $d^2$ over powers of 2, if you want to see what is really possible. I'm totally confused about what you are reporting. Can you enlighten me what you are trying to do, really? | |
| Mar 23, 2023 at 23:14 | comment | added | Yaroslav Bulatov | So the error looks pretty good, eval here. One thing I realized when looking at concrete formula for sum of entries of (1+z A)^{-1} is that it has the form of Heaviside formula when $A$ is diagonal + rank1, so a more direct method for exp(-z A) appears possible | |
| Mar 23, 2023 at 2:39 | comment | added | fedja | @YaroslavBulatov OK, keep me updated (especially if you face some unexpected issues or inexplicable aberrations: the technique is theoretically sound but we both know that implementing something in a stable and reliable way on a real machine has its own challenges :lol:) | |
| Mar 23, 2023 at 2:18 | comment | added | Yaroslav Bulatov | Didn't forget about it, a naive port worked and the underlying math is really amazing so catching up more on background. Plan to clean up and do some evals tomorrow. This method appears to extend beyond DPR1 matrices, any matrix with cheap to compute resolvent appears to benefit from this trick! | |
| Mar 22, 2023 at 11:25 | comment | added | fedja | @YaroslavBulatov Let me know how your implementation worked. I'd like to add a few more final comments about this problem, but I'd rather not bother you with them until you finish the implementation and test it to your satisfaction or otherwise :-) | |
| Mar 21, 2023 at 19:51 | comment | added | fedja | "you are measuring it for scalar arguments and saying it transfers to generalized arguments" Not sure what exactly you mean by that, but if it is "the relative error of the sum of positive terms does not exceed the maximum of the relative errors of individual terms", then yes. The formula for $E$ gives you the worst case scenario (i.e., if the truth is $T$, you computed estimate will be between $T$ and $(1+E)T$), so I would write a program that allows you to choose the degree $n$ of the polynomial (that will require precomputing and storing several sets of complex roots and derivatives). | |
| Mar 21, 2023 at 19:45 | comment | added | Yaroslav Bulatov | For the forward I mean suppose we are able to evaluate exp(x) cheaply , but need (1+x)^{-1}. For the precision, you are measuring it for scalar arguments and saying it transfers to generalized arguments, right? | |
| Mar 21, 2023 at 19:30 | comment | added | fedja | @YaroslavBulatov Great. Just be careful with precision issues! And what exactly do you mean by "an analogous trick of the forward Laplace?" (I.e., what should it accomplish?). Also, have you understood the theoretical relative error estimates I provided? | |
| Mar 21, 2023 at 19:21 | comment | added | Yaroslav Bulatov | Ok I understand it enough now to reimplement from scratch.and do some evals, pretty neat trick. I understand that you are doing something in the spirit of Heaviside formula - constructing a cheap rational expression to write nverse Laplace transform in terms of few samples. I'm curious, is there an analogous sampling interpolation trick for the forward Laplace? | |
| Mar 21, 2023 at 6:22 | comment | added | fedja | @YaroslavBulatov one over the polynomial. The bottom half-plane roots are complex conjugate, so if $A$ is the sum over the top half plane roots, then the sum over the bottom half plane roots is $\bar A$, and instead of honestly computing it again and writing the full sum as $A+\bar A$, you can just take twice the real part, which reduces the running time twice. Note that I used even $n$, so there is no (negative) root on the real line (otherwise it would have to be treated separately). The trick is just partial fraction decomposition. Recall those boring rational function integration drills. | |
| Mar 21, 2023 at 6:13 | comment | added | Yaroslav Bulatov | Thanks for the notes, starting to understand now. So your representation rewrites polynomial in terms of terms that look like resolvents. Does this trick have a name? I don't understand why you only need roots in top half plane | |
| Mar 21, 2023 at 3:27 | comment | added | fedja | @YaroslavBulatov There are some introductory lectures online like www2.math.upenn.edu/~kirillov/MATH548-F07/Lect1.pdf which may help though the best way to understand it is just to try to check the corresponding statements yourself :-) | |
| Mar 21, 2023 at 3:14 | comment | added | fedja | @YaroslavBulatov Have you seen the polynomial operator calculus? If yes, just multiply by the common denominator. Or, for self-adjoint operators, just diagonalize and check everything for the diagonal case. | |
| Mar 21, 2023 at 3:05 | comment | added | Yaroslav Bulatov | Haven't seen this rational calculus trick before, looking for some intro materials | |
| Mar 21, 2023 at 2:05 | comment | added | fedja | @YaroslavBulatov It comes straight from the previous displayed formula: once we have an identity for the rational functions, we can apply it to the operators (rational calculus is always available with one over becoming the inverse of the corresponding operator). Note again that there is a misprint on the RHS (see my comment above) | |
| Mar 21, 2023 at 1:51 | comment | added | Yaroslav Bulatov | I see thanks. I'm trying to understand how you are connecting resolvent to matrix exponential. I'm guessing it's the equation right after "Thus we get"....where does this equation come from? | |
| Mar 21, 2023 at 1:29 | comment | added | fedja | @YaroslavBulatov $\alpha$ and $s_q$ are not related to $P$, but to just finding $w(z)=(I+zV)^{-1}w$ for a fixed $z$). $P$ and $F$ appear later and use this computation to find $w(-r_j^{-1}t/2)$ (another idiotic misprint: $t/2$ escaped from the argument of $w$) in the formula for $P(tV/2)^{-1}w$. | |
| Mar 21, 2023 at 1:09 | comment | added | Yaroslav Bulatov | I'm a little confused about how you are using $\alpha$ and $s_q$ It seems to be the only place that uses $h_i$, but I don't immediately see the connection to $P$. All other mentions of $P$ and $F$ in your description use eigenvalues/eigenvectors of $V$, but those are not efficiently computable, so trying to figure out what $s_q$ are for | |
| Mar 21, 2023 at 0:41 | comment | added | fedja | @YaroslavBulatov Also there are still a few stupid misprints like $F(w)$ instead of the correct $F(t)$, but I do not want to bump the thread to the front page again now, so I hope they are not too confusing :-) | |
| Mar 21, 2023 at 0:13 | comment | added | fedja | @YaroslavBulatov Well, it is pretty much equivalent to my formulas, so I don't think you'll get away with one summation before finding the entries of $w(z)$. Also, if we take into account that addressing the element of an array is like 8 multiplication (or, at least, that is how it was when I learned programming), I guess that won't help anyway. If the code is "faster than you thought possible", try to increase the degree of $P(t)$ a bit. Note that it doesn't use the special structure of the data except the positivity of the entries of $h$. What is the largest $d$ you really want to handle? | |
| Mar 20, 2023 at 23:52 | comment | added | Yaroslav Bulatov | still parsing through the details, but I can run the code and it's faster than I thought was possible. BTW, there's a simple expression for resolvent of $V$ in terms of resolvent of $A$ and resolvent of $h\otimes h$, I wonder if that could be useful here | |
| Mar 20, 2023 at 0:08 | history | edited | fedja | CC BY-SA 4.0 |
edited body
|
| Mar 19, 2023 at 19:23 | history | edited | fedja | CC BY-SA 4.0 |
added 4 characters in body
|
| Mar 19, 2023 at 19:02 | history | edited | fedja | CC BY-SA 4.0 |
added 1 character in body
|
| Mar 19, 2023 at 14:57 | comment | added | fedja | @YaroslavBulatov Here is the underlying "mini-theory". I hope it gives you enough speed in your computations. Also make sure that you have enough precision: I wouldn't dare to raise a number to the power $d^2=10^8$ with just regular double type. Let me know how it works out and feel free to ask as many questions as you want :-) | |
| Mar 19, 2023 at 14:53 | history | edited | fedja | CC BY-SA 4.0 |
added 7477 characters in body
|
| Mar 19, 2023 at 6:52 | comment | added | Yaroslav Bulatov | Thanks for doing this work, I'll dive deeper into this on Monday. The concrete $h$ I'll try to evaluate it on is for $h=1^{-p},2^{-p},\ldots,d^{-p}$ with $p=1.1$ and $d=10000$. Looking forward to your explanation! | |
| Mar 19, 2023 at 6:11 | history | answered | fedja | CC BY-SA 4.0 |