Why does fplot think my function is not vectorized

Executing the following as a script in R2018a

a=ones(1,10);
b=a;
b0=1;

fun=@(x) myFourier(x,a,b,b0);
fplot(fun);

function f=myFourier(x,a,b,b0)

   n=numel(a);
   arg=x(:).*(1:n);
   
   f=b0+sin(arg)*a(:)+cos(arg)*b(:);
   f=reshape(f,size(x));
   
end

results successfully in a plot, but throws warnings (EDIT: and results in non-vectorized execution!!!)

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your
function to return an output with the same size and shape as the input arguments. 
> In matlab.graphics.function.FunctionLine>getFunction
  In matlab.graphics.function.FunctionLine/updateFunction
  In matlab.graphics.function.FunctionLine/set.Function_I
  In matlab.graphics.function.FunctionLine/set.Function
  In matlab.graphics.function.FunctionLine
  In fplot>singleFplot (line 234)
  In fplot>@(f)singleFplot(cax,{f},limits,extraOpts,args) (line 193)
  In fplot>vectorizeFplot (line 193)
  In fplot (line 163)
  In test (line 7) 

But the input function is properly vectorized as is easily verified in simple tests,

 >> fun(rand(1,5))
  
  
 ans =

    3.4371    2.5677   14.1698    2.4490    1.0828
    
  
 >> fun(rand(5,1))
 
 ans =

   14.1745
    8.3619
   -0.2579
    1.7024
    1.5748

Why the warnings, then?

Answer 1

[Walter Roberson]

f1 = fun(1);
f2 = fun([1 2]);
f1 - f2(1)
ans =

      4.44089209850063e-16
      

Your code produces different outputs for the same input, which violates the rule that the calculations must be independent.

[Matt J]

As you can see though, it is a difference of 4.4689e-14 percent, easily attributable to floating point noise. And there's nothing either in the warning messages or fplot documentation that says that should be an issue.

[Walter Roberson]

fplot(@(x) x .* (1 + (length(x)>1)*eps))

generates the warning, so regardless of whether it "should", it does check for equality.

[Matt J]

Tech support got back to me and confirmed Walter's explanation of the problem, and said they will look into remedies for future releases. I passed on to them my suggestions for fixes as well.

Answer 2

[dpb]

Because fplot is just looking at the source code, not the output and it doesn't recognize the reshape operation at the end--only that there are no "dot" operators in the evaluation from which it infers (usually correctly, but as you've shown not infallibly) that the function isn't "vectorized".

Whether there's any chance of being able to get this one right without way more parsing logic than would want to use for performance is, I'm guessing, pretty small. You can't just presume that if the user has a reshape call in the function that always works correctly, either.

[Matt J]

If that's true, it seems like an inexplicably complicated and unreliable approach for the code designers to take. Why wouldn't they implement along the lines of,

function h=fplot(fun,bounds)

        x=linspace(bounds(1), bounds(2),1000);
        
        try %assume vectorized
            
           y=fun(x); 
           
           if ~isequal(size(x),size(y))
            error()
           end
            
        catch ME  %not vectorized - evaluate element by element
        
            y=nan(size(x));
            for i=1:numel(x)
                y(i)=fun(x(i));
            end
            
        end
        
        h=plot(x,y,'-');

end

[Matt J]

Also what if fun() were implemented as a MEX file? How would it parse the source code then?

[Walter Roberson]

The code does not look at the source for the function to make this determination.

In older versions the code ran specific tests such as passing in a vector and checking that the size returned was the same. In current versions, unfortunately the tests were moved into a .p file so we do not know how it works.

[Matt J]

The code does not look at the source for the function to make this determination.

The test below supports the notion that fplot does in fact examine the source code (but why????) .

All the isequal() tests return true on my machine (this is in R2018b now), so both versions of myFourier are returning the exact same numerical output for the exact same numerical input. However, only the first version triggers the warnings.

clc

a=ones(1,10);
b=a;
b0=1;
fun1=@(x) myFourier1(x,a,b,b0);
fun2=@(x) myFourier2(x,a,b,b0);

h=fplot(fun1); x1=h.XData;
h=fplot(fun2); x2=h.XData;

isequal(x1,x2)
isequal(fun1(x1), fun2(x1))


function f=myFourier1(x,a,b,b0)
   n=numel(a);
   arg=x(:).*(1:n);
   
   f=b0+sin(arg)*a(:)+cos(arg)*b(:);
   f=reshape(f,size(x));
   
end



function f=myFourier2(x,a,b,b0)
   n=numel(a);
   arg=x(:).*(1:n);
   
   f=b0+sum(sin(arg).*a(:).',2) +sum(cos(arg).*b(:).',2);
   f=reshape(f,size(x));
   
end

[dpb]

I hadn't gone back and looked at implemenation, granted, so I should have noted was a presumption rather than stating as fact...mea culpa.

From behavior and my thinking of how I would implement the warning, I figured "since everything in Matlab is a function" it could look at either the m-file or if compiled the reference to a dot library routine to determine if dot functions were used.

Without something like that, it seems that such cases as the above tests wouldn't trigger the warning is just an observed behavior for which I didn't otherwise have a ready explanation.

[Walter Roberson]

The test below supports the notion that fplot does in fact examine the source code (but why????) .

aa = fun1(1);
bb = fun2(1);
aaa = fun1([1 2 3]);
bbb = fun2([1 2 3]);
aa-bb
aaa-bbb

ans =

4.44089209850063e-16

ans =

0 0 0

fun1 and fun2 produce the same values as each other for input [1 2 3]. However, fun1 and fun2 produce different values from each other for scalar input [1] . The value produced by fun1(1) is different than the first entry in fun1([1 2 3]) but the value produced by fun2(1) is the same as the first entry in fun2([1 2 3]) . Therefore as far as fplot is concerned, fun2 acts consistently and fun1 does not.

[Matt J]

@Walter, that's a good observation, but even aside from the fact that the difference between aa and bb is just floating point noise, it would only explain the warnings I am seeing if fplot performs exactly the test that you've shown with exactly the same numbers. With these alternative numbers, for example, both functions are in agreement:

aa = fun1(10);
bb = fun2(10);
aaa = fun1([10 11 12]);
bbb = fun2([10 11 12]);
isequal(aa,bb) % = 1
isequal(aa,aaa(1))  %=1
isequal(aaa,bbb)  %=1

[dpb]

'Tis interesting puzzle; one for Yair Altman and Undocumented MATLAB, maybe. :)

I'll just comment my presumption was that probably since the parser has to have "compiled" the function anyway, that it is looking at the internal code generated rather than actually reading the source file itself again...but again, that's all just presumption/guessing about how it seems to detect and be sensitive to whether one has used the dot functions or written code to produce the same result without.

[Walter Roberson]

fplot with default range tests with [0] then [1 2 3] then [1] then [2] and then makes a decision about whether vectorization was correct or not.

Probe code attached.

[dpb]

Good sticktoitniveness, Walter! :)

Makes more sense than my supposition altho is a lot of overhead for what seems little gain to end user for most part...since it doesn't quit but just gives a warning whatever cost there is in the multiple calls to the function to built the plot it goes ahead and does, anyway.

If it gave the warning and hints on how to fix and terminated would seem to make more sense from the efficiency standpoint if that's the intent but then that's somewhat rude, too.

[Walter Roberson]

It is common for people to write expressions that do not give correct results when used with vectors. Often the expressions do not permit vectorized input, such as if it uses ^ instead of .^ . But it is also common for users to use / instead of ./ and in that case the expression can give an output that is the correct size but is wrong . Any difference between one-by-one calculations and vector calculations must be assumed to be an error in the algorithm as applied to vectors.

[dpb]

Agreed, but the warning isn't implying the functional value is wrong; only that it could be calculated more efficiently if it didn't have to loop through the expression...so if that's what they're trying to catch they don't do a good job of explaining the possible problem in the resulting error/warning message.

[Matt J]

Worse than that, because fplot incorrectly judges the function to be nonvectorized, it forms the plot by executing the function in non-vectorized fashion forcing the user to endure slow behavior. I think they should let fplot() operate like integral() already does. Let the user specify whether execution will be vectorized or not, and take responsibility for the result.

[dpb]

This could well be the basis for an enhancement/quality of implemenation improvement request.

Answer 3

[Catalytic]

Here's a workaround to trick fplot into doing the right thing.

Answer 4

[Yair Altman]

As far as I can tell (never mind exactly how), internally a check is made whether the output of fun(1:3) is exactly equaln to the output of [fun(1),fun(2),fun(3)]. Even a tiny FP eps difference will cause the warning to be evoked. If you want to disable the error in run-time, run the following command:

warning off MATLAB:fplot:NotVectorized

[Matt J]

Thanks, Yair. But disabling the warning still leaves the bigger problem of execution efficiency. fplot will still go ahead and evaluate fun in non-vectorized (therefore slow) fashion based on the fun(1:3) test, even with the warning disabled.