Funny pictures of sexps
#31
(08/26/2009, 05:46 PM)bo198214 Wrote: Well my thoughts behind that statement was like this:
The matrix power method applied to fixed points is the regular iteration at that fixed point.
Also I would assume the change of function to be continuous in the development point.
Say we have to real fixed points x0 and x1 and going from x0 to x1 as development point of the matrix power sexp. The regular iteration at x0 would slowly transform into the regular iteration at x1, which we know are different (except for fractional linear functions).

As I consider the intuitive iteration as kinda opposite to the matrix power iteration. I would guess (though this is not supported by numerics yet) that the closer to the fixed point I choose my development point x0 the closer is the resulting islog to the regular slog at the fixed point. (To show that the islog to base sqrt(2) is different from the rslog one could consider the half-iterate isexp(1/2+islog(x)) and show that it has a singularity at 2 which is not true for the regular slog). If not so I would wonder how the islog decides which fixed point to take for being regular at.

...

Well if you have different evidence of any point in my chain of conclusion I would like to hear. My understanding of the complex behaviour of the islog does not suffice that I can support your statement of the spectacular different behaviour at the fixed points. Perhaps you can illustrate with some pictures or thought experiments.
The answer can be found by considering the intuitive slog first (because I don't really see how to get to the intuitive slog from the regular slog by reasoning alone, but perhaps you or sheldon can help me here?).

Start with the islog, developed at 0. Create a straight line, from the origin to the upper primary fixed point, and call this line L1. If we calculate the islog of all the points of L1, we will trace out a curve that goes up and slightly to the left, that eventually straightens into a straight line, with slope equal to about -4.2034 for base e. Call this image L2.

If we travel along L2, and perform the isexp function, we should get back an image of the original line, L1.

However, if, at some point on the image L2, we instead travel in a straight line in the positive real direction (to the right), then as we perform the sexp operation (to get a new image L3), you will have something that resembles the regular sexp developed at the fixed point. But in order to get back to "the action", which I consider anything farther than a unit's distance from the fixed point, you must travel very far to the right. As you do so, you wind around the fixed point a great many times, and this puts you well outside the principal branch of the islog. Indeed, the islog in this particular branch will look pretty much exactly like the rslog.

Note that in the limit as you travel up L1 (perfoming the islog to get an image L2, then picking a point on L2 and performing the isexp while going to the right to get L3) all the way to the fixed point, you get exactly the regular sexp, developed at the fixed point. So the regular sexp at the upper fixed point is what lies "beyond" the top of the sexp developed from the inverse of the intuitive slog. The regular sexp at the lower fixed point is what lies "beyond" the bottom of the intuitive sexp.

The reason that the regular sexp and the intuitive sexp seem different is that, if you followed my instructions, you were travelling up and to the left of the origin as you approached the fixed point. When we think of the regular sexp, however, we view it from somewhere that is closer to "the action" (more than one unit distance from the fixed point), so it's developed probably from somewhere up and to the right, so that as we trace out in the imaginary direction, down towards the origin, we get to "the action" well before we get back to the principal branch of the islog. Thus we can't directly see the relation, going from regular to intuitive. But the relation is obvious (to me anyway) in the other direction.

It is easy to "get to" the regular sexp from the intuitive sexp, by this limiting behavior, and hence to the rslog from the islog. But I haven't been able to reverse the process, at least not in any thought experiment. Perhaps you or Sheldon can help me figure out a way to get "back"?
~ Jay Daniel Fox
#32
(08/26/2009, 08:08 PM)bo198214 Wrote: Are you talking about recentering the Abel function or about my method of shift-conjugation of the original function? Computing the Carleman matrix of exp at some other place than 0 is unexpectedly more time consuming, as I already mentioned. 100 terms take perhaps 10-20 min at my place. While 400 at 0 takes perhaps 5min (without your acceleration method of course)
How are you computing the Carleman matrix? I use iterative matrix multiplication, so it takes only slightly more time to move the development point. I've been doing 250 term matrices in only a couple minutes, give or take. Here is some SAGE code I've been using.

Code:
def Carleman(field, expr, x, size):
    carl = matrix(field, size+1, size+1);
    mulm = matrix(field, size+1, size+1);
    cof = taylor(expr, x, 0, size).coeffs();
    for rrow in xrange(size+1):
        for ccol in xrange(len(cof)):
            col2 = cof[ccol][1]+rrow;
            if col2 < size+1:
                mulm[rrow, col2] = field(cof[ccol][0]);
    vec = vector(field, size+1);
    vec[0] = field(1);
    carl[0] = vec;
    for ii in xrange(size):
        vec = vec * mulm;
        carl[ii+1] = vec;
    return carl;

# Small test run
Flt = ComplexField(27);
size = 5;
shft = Flt(1, 0);
x = var('x');
Cf = Carleman(Flt, e^(x+shft)-shft, x, size)
print Cf;
# Calculate the Abel function
Flt = ComplexField(512);
size = 250;
shft = Flt(1, 0);
time Cf = Carleman(Flt, e^(x+shft)-shft, x, size);
# time CI = Carleman(Flt, x, x, size);
CI = MatrixSpace(RR, size+1, size+1).identity_matrix();
M = Cf.submatrix(1,0,size,size) - CI.submatrix(1,0,size,size);
time A = vector(Flt, M.solve_left(vector(Flt, [1] + [0 for kk in xrange(size-1)]))[0]);
print A.change_ring(CC)[0:20];

Output:
Code:
[  1.000000          0          0          0          0          0]
[  1.718282   2.718282   1.359141  0.4530470  0.1132617 0.02265235]
[  2.952492   9.341549   12.05983   8.945981   4.699514   1.925110]
[  5.073214   24.07712   50.12800   62.18783   53.35036   34.82992]
[  8.717212   55.16170   158.4776   278.1409   340.3287   314.9432]
[  14.97863   118.4792   434.1022   987.6331   1582.217   1922.546]
Time: CPU 52.15 s, Wall: 52.51 s
Time: CPU 20.40 s, Wall: 20.41 s
(0.915948856840961, -0.208605479751221, -0.0545257718234947,
0.0713217792410870, -0.0199817585842940, -0.0110084938578771,
0.0119893515581789, -0.00267961632960444, -0.00264867908203897,
0.00236244492078691, -0.000354279433811006, -0.000647870949133354,
0.000488726691175665, -0.0000401511878585745, -0.000157929045310358,
0.000108539033152715, -2.24888735429443e-6, -0.0000505453252195363,
0.0000483511387610123, -0.0000269817432526313)

Note that I used a ComplexField. A RealField would have sufficed and would have been faster, but this shows that we can now easily step away from the real line if we wanted to.
~ Jay Daniel Fox
#33
(08/26/2009, 08:08 PM)bo198214 Wrote: As far as I remember you only suggested that for recentering. Which I dont do, I thought that became clear now. The development slog at 0, sexp at -1, behaves like expected with convergence radius 1. So there is no "of course" if I do the same procedure just with a conjugated function.
Yes, sorry, I sometimes am not clear, or perhaps ambiguous is the correct term. There is another reason you must truncate, regardless of whether you've recentered, which is that only the first so many terms of the calculated slog taylor series will be "accurate". As such, you can only use that many terms when reverting the series. As I discussed when originally analyzing Andrew's slog, only about the first half of the terms will be accurate.

To be more precise, for the unaccelerated islog, only about the first half of the terms will seem to converge on the correct values, while the rest will be quite wrong, though of roughly the correct magnitude. (To see this, calculate a system twice as large, and note that only the first half of the original terms will roughly match the first quarter of the new terms.)

For the accelerated slog, the values will of course converge on the taylor series for logarithmic singularities at the fixed points, but only about the first half of the terms of the "residue" (as I called it) will seem to converge on the correct values. Again, to see this, calculate a system twice as large, and note that the entire original series will seem to converge on the first half of the new series, but when you look at the "residue", only the first half of the original series will roughly equal the first quarter of the new series.

(And in all cases, in the absence of a proof that the infinite system is solvable, "correct" means the converging values that we seem to get as we use larger and larger matrices, so we can only analyze matrices much smaller than the largest we have yet managed to solve.)

For the unaccelerated islog, I haven't analyzed how many terms you can successfully revert, though it's probably much less than half. For the accelerated islog, it's just about half, give or take, and I don't have my original notes on it, so I can't be more specific without redoing my calculations. It may have been more than halfway and I was just being conservative, but I remember for sure that 900 terms out of 1200 just got way off, and that was with the accelerated solution.

Quote:
(08/26/2009, 07:21 PM)jaydfox Wrote: Count "numerically" as a singularity? It's not a singularity, so at best it means that you just need more terms in the power series.
Of course that was what I was saying. Little changes in the intput cause big changes in the output that means increasing the precision, doesnt it?
Ah, sorry.
~ Jay Daniel Fox
#34
(08/26/2009, 08:38 PM)jaydfox Wrote: How are you computing the Carleman matrix?
Actually I use a direct formula, I thought that was the fasted method. But it seems that I am wrong:
Code:
C = Matrix([ [log(b)**n/factorial(n)*sum([binomial(m,k)*k**n*(b**x0)**k*(-x0)**(m-k) for k in range(m+1)]) for n in range(N)] for m in range(N)])
(Where N is the size, b is the base and x0 is the shift)
I also always used precision 1000 but that does not stretch the time in that amount.
I really dont know why it takes so long.

Why is in your Carleman matrix output only little precision while you put in 512?

Quote:I use iterative matrix multiplication,

I am not sure how it works. Can you give just the basic idea?
#35
(08/26/2009, 09:15 PM)bo198214 Wrote: Why is in your Carleman matrix output only little precision while you put in 512?

Oh I see because ChangeRing
#36
(08/26/2009, 09:15 PM)bo198214 Wrote: I am not sure how it works. Can you give just the basic idea?

Ok. I think I also understood that now, the lines are iteratively computed and perhaps that saves lots of time because the result from the previous step is used.
#37
(08/26/2009, 09:15 PM)bo198214 Wrote: I also always used precision 1000 but that does not stretch the time in that amount.
I really dont know why it takes so long.
Well, as far as the amount of precision, this needs to be tuned to the problem being investigated. For the slog, I found that I need about 1.4 times as many bits as terms (which oddly enough is close to the magnitude of the radius of convergence, but this might be a coincidence), just to get numerical stability, and then however many additional bits of accuracy I want. So if I want 200 bits of accuracy, and I'm using 250 terms, I need about 550 bits, maybe a bit less (520?), so I just used 512 to pick a nice multiple of 64.

I'm not sure if the same 1.4x rule applies for your recentered function; it could be 1.3 or 1.5 or something like that, so a bit of empirical investigation is required before you can safely tune to tight bounds. (Yes, I know your method is not recentered in the sense of recentering a polynomial, but it's still recentered in terms of what the output is compared to the input).

Quote:
Quote:I use iterative matrix multiplication,

I am not sure how it works. Can you give just the basic idea?
Well, the Carleman matrix consists of rows that are powers of the Taylor series in question (iterated multiplication, not iterated composition). To calculate the multiplication, we need to perform "long multiplication", which can be done with a matrix (the matrix called "mulm" in my code).

Probably the best way to see this is with a simple code example:

Create a matrix:
Code:
var('a,b,c,d,f,g')
M1 = [[a, 0, 0],
     [b, a, 0],
     [c, b, a]]
M2 = [d, f, g]
M3 = [];
# Multiply the matrix by the vector (can't make symbolic matrices
# in SAGE, so I'm faking it here with lists)
for kk in xrange(3):
    M3.append(sum(M1[kk][jj] * M2[jj] for jj in xrange(3)));
print M3
f1 = a + b*x + c*x^2;
f2 = d + f*x + g*x^2;
f3 = (f1*f2).expand().series(x, 3);
print f3;

Output:
Code:
[a*d, a*f + b*d, a*g + b*f + c*d]
(a*d) + (a*f + b*d)*x + (a*g + b*f + c*d)*x^2 + Order(x^3)
~ Jay Daniel Fox
#38
(08/26/2009, 09:55 PM)jaydfox Wrote: Well, the Carleman matrix consists of rows that are powers of the Taylor series in question (iterated multiplication, not iterated composition).

The matrix application seems a bit obfuscated to me. Is there any reason why you dont do:
Code:
def psmul(A,B):
    N = len(B)
    return [sum([A[k]*B[n-k] for k in xrange(n+1)]) for n in xrange(N)]

C = Matrix(R,N)
row = vector(R,[1]+(N-1)*[0])
C[0] = row
for m in xrange(1,N):
    row = psmul(row,coeffs)
    C[m] = row
where coeffs is the sequence of coefficients and R the corresponding Real/ComplexField?

Btw. the code is now fast. Thanks Jay.
#39
(08/26/2009, 10:26 PM)bo198214 Wrote: The matrix application seems a bit obfuscated to me. Is there any reason why you dont do...
Well, I've run into problems with creating matrices and vectors from lists, especially long lists. I've seen code take minutes with lists that takes only seconds if explicitly set with for loops. The problem only pops up in certain circumstances, but I've been bitten by performance issues enough times that I just start making everything explicit, just to be sure. There are a few circumstances where the implicit method is faster, but I usually start with implicit code to make sure my methods are sound, then explicit code for speed, and then if I have time, I try to see where I can sneak some implicit code back in for even more speed (or memory savings).

Also, the explicit code is definitely more memory-friendly in most cases, and when I try to push these methods to their limits, memory-friendly code is a must.

Update:
(08/26/2009, 10:26 PM)bo198214 Wrote: Btw. the code is now fast. Thanks Jay.

Yes, your version is even faster, more than twice as fast as mine (half as many multiply-adds during the power series multiplication)!
~ Jay Daniel Fox
#40
(08/26/2009, 08:11 PM)jaydfox Wrote: The answer can be found by considering the intuitive slog first (because I don't really see how to get to the intuitive slog from the regular slog by reasoning alone, but perhaps you or sheldon can help me here?).

Start with the islog, developed at 0. Create a straight line, from the origin to the upper primary fixed point, and call this line L1. If we calculate the islog of all the points of L1, we will trace out a curve that goes up and slightly to the left, that eventually straightens into a straight line, with slope equal to about -4.2034 for base e. Call this image L2.

....

I dont think one can not understand this without having taken a proper look at the pictures.
Fortunately what you (also previously) describe does not only apply to the intuitive slog but I remember the pictures that Dmitrii made about his Cauchy slog.
(They are scattered over the forum. Particularly this image depicts I think L2, as the left border. It depicts cslog(G) where G is the region bounded by L1 and exp(L1), though L1 extends here also to the lower fixed point. This region is particularly also used by Kneser in his construction which though is not of computational nature.)

So remembering these pictures I can follow you at least a bit.
If I believe what you write I am asking myself where the specific *intuitive* slog is relevant. I guess you could do that construction with every slog and always get as a limit the regular slog. Because you would continue every slog that is defined on G to the whole plane except to the primary fixed points and cuts and somehow it seems as these coarse properties suffice for your construction instead of it relying on the fine grained structure directly on G.

But I admit I am not so deep into that construction (particularly I didnt optically verify your described resemblence to the rsexp) so you tell me!


Possibly Related Threads…
Thread Author Replies Views Last Post
  Roots of z^z^z+1 (pictures in MSE) [update 8'2022] Gottfried 5 8,789 08/30/2022, 02:08 AM
Last Post: JmsNxn
  On n-periodic points of the exp() - A discussion with pictures and methods Gottfried 1 4,673 06/10/2020, 09:34 AM
Last Post: Gottfried



Users browsing this thread: 3 Guest(s)