I'm trying to construct a set of orthogonal polynomials, starting from a specially-prepared initial polynomials:

nMax = 3;
initPolys = Table[1 + R/2 x - R/2 x^2 - x^n, {n, 2, nMax}];
orthoPolys = Orthogonalize[initPolys, Simplify @ Integrate[#1 #2, {x, 0, 1}] &]

The result of Orthogonalize grows rapidly in size as nMax grows to as small values as 7. Yet, if I use Simplify on the output, the expressions appear quite "compressible".

The problem is that this simplification only happens after the whole set has been generated, and in the Gram-Schmidt method used by Orthogonalize, the large unsimplified expressions are being used, thus leading to unacceptable slowdown. E.g. nMax=11 has already been running for 10 hours without any result (so I aborted it).

Thus my question is, how can I make Orthogonalize simplify each new vector before using it in further calculations? Or do I have to make my own implementation of Gram-Schmidt process to do this?

  • 1
    $\begingroup$ Oddly enough, if you route your answer back through the Orthogonalize (so orthogonalize the already orthonormal set, Orthogonalize[orthoPolys, Simplify@Integrate[#1 #2, {x, 0, 1}] &]) it takes longer than the original Orthogonalize. $\endgroup$ – MikeY Jan 20 '19 at 23:46
  • $\begingroup$ @MikeY that's expectable: the second time it has to integrate&simplify larger expressions when calculating dot products. $\endgroup$ – Ruslan Jan 21 '19 at 5:53
  • $\begingroup$ For a bit, I was thinking that you could orthogonalize the smaller problem then just build up repeatedly from there. Might still be a solution based on that approach, if you can “tell” the algorithm the earlier terms are already orthogonal. $\endgroup$ – MikeY Jan 21 '19 at 13:28

Recalling that the coefficients of the product of two polynomials is the same as the convolution of the coefficients of the factors, one can easily write a more efficient version of the required inner product:

iprod[f_, g_, x_] := With[{c = Simplify[ListConvolve[CoefficientList[f, x], 
                                                     CoefficientList[g, x], {1, -1}, 0]]}, 

Evaluating Orthogonalize[Table[1 + R/2 x - R/2 x^2 - x^n, {n, 2, 9}], iprod[#1, #2, x] &] is now much faster, tho the resulting polynomials need to be processed further with (Full)Simplify[].


I did some timing experiments with your example problem. Here is my code:

N[Round[#1[[1]], 10^-3], 3] & /@ Table[ initPolys =
  Table[1 + R/2 x - R/2 x^2 - x^n, {n, 2, nMax}]; Timing[
  orthoPolys = Simplify@Orthogonalize[initPolys, 
  Simplify@Integrate[#1 #2, {x, 0, 1}] &]], {nMax, 9}]

I ran this code on an old PC with some older versions of Mathematica:

 6.0: {0, 0.0320, 0.0460, 0.125, 0.281, 0.499, 0.780, 1.14, 1.75}
 7.0: {0, 0.0160, 0.0620, 0.141, 0.280, 0.500, 0.842, 1.26, 1.84}
 8.0: {0, 0.0160, 0.1250, 0.577, 0.936, 2.31, 6.16, 8.81, 15.8}
 9.0: {0, 0.125, 0.967, 2.40, 7.33, 14.6, 25.6}
10.2: {0, 0.109, 0.250, 2.42, 8.27, 18.0, 31.4}
11.3: {0, 0.0240, 0.147,0.823, 3.17, 13.3, 29.9}

For versions 9+ I only went up to nMax=7 and the version 11.3 ran on Wolfram Development Platform in the cloud. The actual runtimes are variable depending on the CPU, and so on, but the trend is clear. For versions up to 9 the growth seems to be O(n^3.5) up to O(n^4). For versions 10+ the rate goes up to O(n^5) or O(n^6). This is the reason for your excessive runtimes. An impractical suggestion is to use version 6.0 instead. Another suggestion is to use numerical values for R.

  • $\begingroup$ Numerical values for R will make the resulting polynomials lose orthogonality as we increase size of the basis set. Gram—Schmidt orthogonalization isn't a stable process. Even if Mathematica tried to use a stabilized version of it, this would only delay loss of orthogonality, but not eliminate it. $\endgroup$ – Ruslan Mar 12 '19 at 7:27
  • $\begingroup$ @Ruslan, there are in fact stabilizations of the Gram-Schmidt process that might be more suitable, but they are not yet implemented in Mathematica, and thus need to be programmed. $\endgroup$ – J. M.'s ennui Mar 12 '19 at 7:41
  • $\begingroup$ @J.M.isslightlypensive what stabilizations could be better suited for polynomials than modified Gram–Schmidt? I can't seem to find any at all — as if MGS is the only one known. $\endgroup$ – Ruslan Mar 12 '19 at 9:27
  • $\begingroup$ @Ruslan, I haven't yet found time to read the papers and do experiments, but it appears there is some research in reorganizing the steps in MGS for orthogonalizing polynomials, Moving from GS-type methods, I have seen algorithms based on "quasimatrices", where the equivalent of a Householder transformation on matrices is done to a set of basis polynomials instead. The orthogonal polynomial literature is really rich! $\endgroup$ – J. M.'s ennui Mar 12 '19 at 9:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.