Vectorizing a function looking for eigenvalues of a matrix

In short, I have a function that defines a matrix based on an input value, and then returns the eigenvalue of that matrix that is closest to 1. I want to plot this for different inputs, but I get the following warning:

'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.'

Is it possible to vectorize this function in order to plot it properly (and without waiting a very long time)?

My code looks somewhat like this (if necessary I can post the entire code but it's pretty long).

eigenvalue = @eigenvalue_func
fplot(eigenvalue)

function eta = eigenvalue_func(z)
    p = 1:1:100;
    q = p;
    weights = rand(100,1); %in reality I'm using a distribution for p, q and the weights based on Gaussian Quadrature
    [Q, P] = meshgrid(q,p);
    [W, ~] = meshgrid(weights, p);
    v_0 = -30./(P.*Q.*(P.^2-Q.^2).*2.*pi.^2).*(Q.*cos(Q).*sin(P)-P.*cos(P).*sin(Q));
    A = -4.*pi.*v_0.*(Q.^2  - z  ).^(-1).*Q.^2.*W;
    for j = 1:100 %the diagonal has to be defined differently since the above returns NaN/Inf (I take the limit here)
        v_0(j,j) = -30./(2.*pi^2.*P(j,j).^2).*(0.5-(sin(2.*P(j,j)))./(4.*P(j,j)));
          A(j,j) = -4.*pi.*v_0(j,j).*(Q(j,j).^2  - z).^(-1).*Q(j,j).^2.*W(j,j);
    end
    eta_full = eigs(A, 20);
    [~, index_eigenvalue] = min(abs(eigs(A, 20)-1));
    eta = eta_full(index_eigenvalue);
end

[Torsten]

At least you should show us which call to the function "eigenvalue_func" generated the error message you mention.

[Bruno Luong]

Torsen it is called internally by

fplot(eigenvalue)

[J. Alex Lee]

I don't think you can: z in the context inside the function has to be a scalar, right? The error would come from fplot wanting to feed it vector z of arbitrary length, so it is not generally compatible to subtract from array of size(Q).

[Ferdy Ellen]

z is indeed a scalar. I was hoping there was some clever way of doing it but all I could think of would be to use a for loop to take the difference elementwise (for z as a vector) but I don't think that will speed things up at all

[Bruno Luong]

The bottleneg is probably calling EIGS (have you made any profiler of your code?) which cannot be vectorized at the moment, and possibly fplot waste a lot of time to find the interval and resolution to plot the function.

If you want faster plot, change fplot and eigs.

For instant it is totally iefficient to call twice eigs with 20 eigenvalues

eta_full = eigs(A, 20);
[~, index_eigenvalue] = min(abs(eigs(A, 20)-1));
eta = eta_full(index_eigenvalue)

where you can get the same result with

eta = eigs(A,1,1); 

[Walter Roberson]

Since eigs is not internally vectorized, it might make sense to define a function that accepted a vector of z values, and used background pool to evaluate multiple z in parallel.

Answer 1

[Walter Roberson]

wrapper = @(z) arrayfun(eigenvalue, z)
fplot(wrapper)