Meaning of 'degree' is erroneous or wrong

[hwborchers]

I did one of my favorite tests, approximating the Runge function with a polynomial of degree 10 (i.e., 11 coefficients) on the interval [-1, 1].
See my discussion on StackExchange "Computational Science" from 2012:
https://scicomp.stackexchange.com/questions/1531/the-remez-algorithm

Doing this with degree = 10, I received the following error message:

    library(minimaxApprox)
    lb <- -1.0; ub <- 1.0
    fnRunge <- function (x) 1/(1 + (5 * x)^2)

    sol <- minimaxApprox(fn = fnRunge, lb, ub, degree = 10)
    ## Error in solve.default(polyMat(x, y, relErr), y) : 
    ##  system is computationally singular: 
    ##  reciprocal condition number = 2.00539e-19

I wondered if I misunderstood the 'degree' parameter, but setting degree = 11 returned a vector of 12 elements, i.e., a polynomial of degree 11.

    sol <- minimaxApprox(fn = fnRunge, lb, ub, degree = 11)
    sol$a
     [1]  9.340771e-01 -1.855657e-15 -1.155302e+01  3.126176e-13  5.917189e+01
     [6] -1.939782e-12 -1.341553e+02  4.249187e-12  1.357960e+02 -3.850310e-12
    [11] -5.022113e+01  1.230144e-12

The highest coefficient is almost 0.0, so forgetting it, we get a vector of length 11, representing a polynomial of degree 10.

This corresponds reasonably well with the result we get when we apply the Remez algorithm in Chebfun in MATLAB or the Remez.jl module in Julia:

    p <- c(0.934077073, 0.0, -11.553015692, 0.0, 59.171892231,
           0.0, -134.155250367, 0.0, 135.795965068, 0.0, -50.221129702)

and the estimated maximal distance is 0.06592293 in all three calculations.

CONCLUSION

I strongly believe you should reconsider the meaning of 'degree' in the call to the minimaxApprox function. It does not coincide with the description "Either a single value representing the requested degree for polynomial approximation ..." on the help page.

And the error for degree = 10 is not correct - there is a solution, only that the highest coefficient is zero.

[aadler]

Hello, @hwborchers. Thank you for the note. Yes, you are correct that sometimes the algorithm will converge when requesting one extra term and that term is almost 0, where it will fail otherwise. I note this in the error message for rational approximation when there is a 0 in the range for the denominator polynomial (https://github.com/aadler/minimaxApprox/blob/master/R/RemezRational.R lines 70–74).

As for the degree 11 vs. 10, I did find in my research/experimentation that when "zapping" small values, the algorithm failed more frequently. So I am loath to simply ignore those terms. But I do know that certain classes of functions should have alternating terms as 0. To build that logic into the algorithm will take more research, study, and understanding on my part. If you have any suggestions as to how to proceed, please let me know.

Most importantly, to address your conclusion, I would say that I use "degree" to mean the maximum exponent of the requested polynomial, which is why it calls for n + 2 basis points. I believe he error is technically valid in that it isn't saying that the requested n-degree polynomial (10 in your case) does not exist, but that the algorithm as implemented failed to find one due to the singular nature of the matrix and the limitations of floating-point implementations. Perhaps a better error message would help? If you disagree, what would you suggest instead? Also, any suggestions you may have about the package in general would be greatly appreciated. Thank you!!

[aadler]

Addressed in commit cbe39d9. Will be part of next release.