recurrence relation for any given 'n'.

How to compute A_j^(2n) for any 'n' using A_j^(2).

Here, n=2,3,4,...; and j=1:n-1.

Any kind of help is highly appreciated.

Thank you.

Answer 1

[Walter Roberson]

If you have the Symbolic toolbox, you could use evalin(symengine) or feval(symengine) to call into MuPAD in order to take advantage of solve() of a rec() recurrence structure

[Walter Roberson]

Guessing that there is a vector of initial constants, A(J)^1, J = 1..n-1, and calling that vector A_1:

function result = Aj_2(j, A_1)
  result = A_1(1:j) * A_1(j:-1:1).';
end

function result = Aj_2n(j, two_n, A_1)
  if two_n == 2
    result = Aj_2(j, A_1);
  else
    result = 0;
    for mu = 1 : j
      result = result + Aj_2n(mu, two_n - 2, A_1) .* Aj_2(j - mu + 1, A_1);
    end
  end
end

You can invoke this by creating a specific vector of values:

function result = A_driver
  N = 7;                  %for example
  A_1 = randi(5, 1, N);   %for example

    result = zeros(1, N-1);
    for j = 1 : N-1
       result(j) = Aj_2n(j, 2*N, A_1);
    end
  end

With large enough values of N, the cost of recursion might add up significantly (though you would probably risk exceeding 2^53, the point at which you stop being able to distinguish adjacent integers.) I was going to show a way in R2017a and later to improve that speed, but the code above is fast enough for the values I tried that it is not worth going to the trouble.

[Venkata]

Dear Roberson,

Thank you for your prompt response. However, I got large values from your computation.

The multiplication in the first function, 'result' should be a vector multiplication not the matrix multiplication. I corrected that, then I received the following values, for N=5 and with A_1 having the actual values.

641584391529617.250000
12978244794078948.000000
147984150516881280.000000
1267037968968080400.000000

However, what I am supposed to get are the following values. These I got from my computation. I will cross check these too.

         915.062500
         3702.055500
        12257.898481
        41869.900040

Here are the A_j^(2) [it is the variable 'result' in your code] values I have used:

           30.250000
           61.191000
          140.720000
          407.410000
         1419.200000

Thank you so much.

[Walter Roberson]

"The multiplication in the first function, 'result' should be a vector multiplication not the matrix multiplication."

No, that is not correct. Notice that I used

A_1(1:j) * A_1(j:-1:1).'

with the second item transposed to create a column vector. That makes it a row vector matrix-multiplied by a column vector, which will have the effect of multiplying the corresponding locations and then summing the result. The code is equivalent to

sum( A_1(1:j) .* A_1(j:-1:1) )

which could also be written

dot( A_1(1:j), A_1(j:-1:1) )

Of these, the form that I used is the fastest, with the explicit sum of .* the second fastest, and dot() about 1/3 slower than the form I used.

[Walter Roberson]

I found that with larger N the recursion was a big problem. Here are versions that have considerably improved performance for larger N:

[Venkata]

Dear Roberson,

Yes, you are right. That was not a matrix multiplication. I made it a column vector by using A_1(:) and that's why got a matrix. My bad. I am sorry.

Nevertheless, I get, 'memoize' is an undefined function or variable error. I am using Matlab-2016a version, but I saw this feature was explained by Loren in 2006 in Mathworks page. So, should I activate some feature to use this facility in my version?

Thank you.

[Walter Roberson]

memoize() was new in R2017a. You can acheive the same effect using tools such as the one mentioned at http://blogs.mathworks.com/pick/2015/03/27/cache-your-nuts-to-eat-later/, but that does say something about writing to disk.

With your N being around 5, the memoization / caching is not important and the original would be fine. I just happened to randomly toss in 17 for N and the result was so slow that I was inspired to try out the R2017a memoization feature. If you are only ever going to be using small values for N then it is not worth the time to create a version with caching for use with earlier MATLAB releases.

[Venkata]

Dear Roberson,

Thanks a ton. That was perfect.

Thank you so much.

For the testing purpose I chose N=5 else there is no limit to it. I want up to N=50 at least.

Now I know how to do it.

I can't thank you less.

Thank you.

[Walter Roberson]

memoization / caching would make a huge difference for N = 50.

If you do not use caching, or do not use a large enough cache, you would probably be better off building upwards, starting from the smaller values and calculating up to the larger values using the already-calculated values.

[Venkata]

Dear Roberson,

Is there any work around for memoization in 2016a?

Seems that I need to run for N=1500. I ran the code for N=500 yesterday evening and even after 12 hrs it's still running.

Please let me know how to achieve memoization in 2016a.

Thanks a lot.

[Walter Roberson]

I tossed together the attached. I did not test it with functions that return multiple outputs.

Aj_2n_better and Aj_2 are the same as last time. memoizefun is new (I just wrote it) and A_driver_cache is almost the same as the previous A_driver_better except for calling memoizefun instead of R2017a's memoize()

Note: you should test that the results are the same as the previous... they should be, but I did not test that.

[Venkata]

Dear Roberson,

Sure, I'll check them.

I owe you big time, Mr. Roberson.

Thanks a ton.