Force imaginary part of complex answer to be postive (efficiently)

I am using fsolve to find the complex roots over a very large parameter space (which it is doing quite well). However, I would like the system to only return the positive pair for neatness sake. Which of these would be more efficient to run (yes, I know this is needless optimization, but I have some time while the code runs and thought I would ask)?

if imag(result) < 0
    result = result - 2i*imag(result); % option 1
    % or
    result = real(result) - 1i*imag(result); % option 2
    % or
    result = complex(real(result), -imag(result)); % option 3
    % or ...?
end

If someone who is knowledgeable could explain why each case is faster/slower, that would be really neat to learn!

Thanks!

Answer 1

[Walter Roberson]

format long g
S = 10000;
result = complex(randn(S,1), randn(S,1));

N = 100;
t1 = zeros(N,1); t2 = zeros(N,1); t3 = zeros(N,1); t4 = zeros(N,1);

for K = 1 : N; start = tic; newresult = result - 2i*imag(result); t = toc(start); t1(K) = t; end
for K = 1 : N; start = tic; newresult = real(result) - 1i*imag(result); t = toc(start); t2(K) = t; end
for K = 1 : N; start = tic; newresult = complex(real(result), -imag(result)); t = toc(start); t3(K) = t; end
for K = 1 : N; start = tic; newresult = complex(real(result), abs(imag(result))); t = toc(start); t4(K) = t; end

subset = [t1, t2, t3, t4];
subset = subset(11:end,:);
mean(subset, 1)

ans = 1×4

1.0e+00 * 4.24888888888889e-05 4.35222222222222e-05 2.19888888888889e-05 2.33222222222222e-05

plot(subset)
legend({'option 1', 'option 2', 'option 3', 'unconditional'})

We can see that using your option 3, complex() is notably faster than your first two options.

The 4th option I put in, also using complex() but also using abs(), is an option that does not require the if first. If you were doing a bunch of these, your three options should all properly be timed inside a loop that did an "if" -- that or should be written in terms of logical masks. But my 4th option using complex() and abs() is just slightly slower than your unconditional complex() ... except the timing for your unconditional complex does not include the tests for negative complex part, so it is rather likely that using the 4th option I show here, and no if, would be faster than testing which entries to convert.

[Sanders A.]

Thanks for such an involved answer!

I think I shall go with option 3 going forward because the system is relatively rarely producing results such that the if-statement evaulates to true (which I assume is a faster check than re-creating the complex number regardless of if it's needed or not).

Thanks!

[Walter Roberson]

The cost of doing the test does indeed add up.

In the below, I used three approaches:

• logical indexing
• for loop with individual tests
• unconditional

The first plot shows the logical indexing together with the unconditional approach. The second plot shows the for loop approach.

You can see that the logical indexing approach is significantly faster than any of the for loop approaches. And you can see that unconditional is the fastest approach.

Also I added an approach using conj(); you can see that it is slower than my unconditional abs() approach, but it is faster than any of the other approaches.

format long g
S = 10000;
result = complex(randn(S,1), rand(S,1));
which = rand(S,1) < 0.1;
result(which) = conj(result(which));

N = 100;
t1 = zeros(N,1); t2 = zeros(N,1); t3 = zeros(N,1); t4 = zeros(N,1);
t1L = zeros(N,1); t2L = zeros(N,1); t3L = zeros(N,1); t4L = zeros(N,1);
t5 = zeros(N,1);


for K = 1 : N; start = tic; newresult = option1(result); t = toc(start); t1(K) = t; end
for K = 1 : N; start = tic; newresult = option2(result); t = toc(start); t2(K) = t; end
for K = 1 : N; start = tic; newresult = option3(result); t = toc(start); t3(K) = t; end
for K = 1 : N; start = tic; newresult = option4(result); t = toc(start); t4(K) = t; end
for K = 1 : N; start = tic; newresult = option1L(result); t = toc(start); t1L(K) = t; end
for K = 1 : N; start = tic; newresult = option2L(result); t = toc(start); t2L(K) = t; end
for K = 1 : N; start = tic; newresult = option3L(result); t = toc(start); t3L(K) = t; end
for K = 1 : N; start = tic; newresult = option4L(result); t = toc(start); t4L(K) = t; end
for K = 1 : N; start = tic; newresult = option5(result); t = toc(start); t5(K) = t; end


subset = [t1, t2, t3, t4, t1L, t2L, t3L, t4L, t5];
subset = subset(11:end,:);
mean(subset, 1).'

ans = 9×1

6.82e-05 7.05222222222222e-05 6.61888888888889e-05 4.78444444444444e-05 0.000672855555555556 0.00102005555555556 0.000336577777777778 0.000217533333333333 2.88888888888889e-05

subplot(2,1,1)
plot(subset(:,[1:4, end]))
legend({'option 1 mask', 'option 2 mask', 'option 3 mask', 'conj mask', 'unconditional'});
subplot(2,1,2)
plot(subset(:,5:8))
legend({'option 1 loop', 'option 2 loop', 'option 3 loop', 'conj loop'})

function newresult = option1(result)
    newresult = result;
    mask = imag(result) < 0;
    newresult(mask) = result(mask) - 2i*imag(result(mask));
end

function newresult = option2(result)
    newresult = result;
    mask = imag(result) < 0;
    newresult(mask) = real(result(mask)) - 1i*imag(result(mask));
end

function newresult = option3(result)
    newresult = result;
    mask = imag(result) < 0;
    newresult(mask) = complex(real(result(mask)), -imag(result(mask)));
end


function newresult = option4(result)
    newresult = result;
    mask = imag(result) < 0;
    newresult(mask) = conj(result(mask));
end

function newresult = option1L(result)
    newresult = result;
    for mask = 1 : numel(result)
        if imag(result(mask))
            newresult(mask) = result(mask) - 2i*imag(result(mask));
        end
    end
end

function newresult = option2L(result)
    newresult = result;
    for mask = 1 : numel(result)
        if imag(result(mask))
            newresult(mask) = real(result(mask)) - 1i*imag(result(mask));
        end
    end
end

function newresult = option3L(result)
    newresult = result;
    for mask = 1 : numel(result)
        if imag(result(mask))
            newresult(mask) = complex(real(result(mask)), -imag(result(mask)));
        end
    end
end

function newresult = option4L(result)
    newresult = result;
    for mask = 1 : numel(result)
        if imag(result(mask))
            newresult(mask) = conj(result(mask));
        end
    end
end


function newresult = option5(result)
    newresult = complex(real(result), abs(imag(result)));
end