Your comments
Oh, I see the difference now, thank you (never check mails on the go...).
So, TWM added conjugation to cross function in R2019b.
In Update 3 they reverted it back to pre-R2019b state (without conjugation).
Will add workaround in next version of toolbox.
Thank you for the report.
Strange, the "cross.m" was not changed in Update 3. What version of toolbox do you use on Windows?
Please update it if it is too old.
Looks like TWM just fixed the bug they introduced in initial release of R2019b.
While I am installing it, could you please share code of cross.m distributed with Update 3 (email or here)?
Any follow up on this?
Thank you for your response.
I have done some quick tests using ill-conditioned problems from NLEVP library. The 'polyeig' provided in our toolbox works quite well on them.
I wrote simple function 'test_polyeig' which returns absolute error of eigen-decomposition for the the problem specified (see its full code at the bottom) for 'double' and extended precision ('mp') cases:
>> test_polyeig('power_plant','double')
ans =
1.8567e+05
>> test_polyeig('power_plant','mp')
ans =
4.091947679918651639316533910037398e-13
>> test_polyeig('cd_player','double')
ans =
1.1410e+04
>> test_polyeig('cd_player','mp')
ans =
4.953565354331632662743256645399015e-15 >> test_polyeig('damped_beam','double') % use 'double' precision ans = 34.2276 >> test_polyeig('damped_beam','mp') % use extended precision ans = 7.547834125835157737664345305865886e-17
Our 'polyeig' gives much more accurate result out-of-the-box. The only possible downside compared to 'quadeig' is that 'polyeig' is (much) slower. But it is perfectly usable for research right now.
Could you please justify why you think "...without 'quadeig' is supported, the toolbox is scarcely useful for a large number of people who work on quadratic eigenvalue problems." Quadratic (and other) eigenproblems can be handled by toolbox right now using 'polyeig'. Or am I missing something?
The code of 'test_polyeig' below. It supports 'double' and 'mp' computations:
function abserror = test_polyeig(name,float_t,varargin)
[COEFFS,~] = nlevp(name,varargin{:});
p = length(COEFFS);
A = cell(p,1);
for k = 1:p
A{k} = cast(COEFFS{k},float_t);
end
n = size(A{1},1);
[X,E] = polyeig(A{:});
error = zeros(n*(p-1),float_t);
for j = 1:n*(p-1)
lambda = E(j);
x = X(:,j);
if lambda ~=Inf
T = A{1};
for k = 2:p
T = T + lambda^(k-1) * A{k};
end
error(j) = norm(T*x,Inf);
else
error(j) = norm(A{1}*x,Inf);
end
end
abserror = norm(error,Inf)
P.S.
To get even higher accuracy, NLEVP code must be converted to support 'mp' computations.
Because right now, NLEVP computes matrices in 'double' precision and then we use 'mp' to compute eigenvalues.
If matrices would be computed with high precision from the onset, than overall accuracy would be higher.
Toolbox uses binary representation of floating-point numbers.
When we extend double precision number with zeroes (in binary representation) - it doesn't refer to zeros in decimal representation (please refer to some guides on floating-point numbers & their formats for deeper understanding).
Do not use mp(1/3) because 1/3 is computed in double precision and then converted to quadruple - final value has double precision accuracy (not the intended quadruple).
The correct way to compute constants to full accuracy is:
>> mp.Digits(34); >> mp('1/3') ans = 0.3333333333333333333333333333333333 >> mp.Digits(50); >> mp('1/3') ans = 0.33333333333333333333333333333333333333333333333333 >> mp.Digits(70); >> mp('1/3') ans = 0.3333333333333333333333333333333333333333333333333333333333333333333333
In this case, double precision constant is not extended by zeros (in binary format) but computed with full accuracy from the onset (within limits of specified precision).
***
Could you please reply to another thread you started yesterday: https://mct.userecho.com/communities/1/topics/179-it-would-be-fantastic-if-quadeig-is-supported
We prioritize only well-motivated requests to add new functionality.
Thank you very much for your request. We have added "quadeig" to our TODO list.
Could you please provide a bit more details regarding your request? This would help us to add support for "quadeig" faster.
- Is there simple example of quadratic eigenvalue problem to demonstrate the need for high-precision (e.g. where double precision fails due ill-conditioning, etc.)?
- In general, what kind of applied & research problems would high-precision "quadeig" enable scientists to solve?
- Are these kind of quad-eig problems common (in research and in applications)?
- What codes researchers in the field are using now? (e.g. Francoise Tisseur's library:https://www.maths.manchester.ac.uk/~ftisseur/misc/)
This information would help us a lot.
Thank you in advance.
Just re-download and use 4.6.4.13348 version.
I have just added support for "cross" from R2019b.
TMW changed the functionality of "cross" function in R2019b for complex vectors.
< R1019b:
% Calculate cross product
c = [a(2).*b(3)-a(3).*b(2);
a(3).*b(1)-a(1).*b(3);
a(1).*b(2)-a(2).*b(1)];
>= R2019b:
% Calculate cross product
c1 = conj(a(2).*b(3)-a(3).*b(2));
c2 = conj(a(3).*b(1)-a(1).*b(3));
c3 = conj(a(1).*b(2)-a(2).*b(1));
if iscolumn(a) && iscolumn(b)
c = [c1; c2; c3];
else
c = [c1, c2, c3];
end
Second one seems to be mathematically more correct. But if we follow this, toolbox would fail on older MATLAB versions.
Have to rewrite the function to support both :(.
In any case, just ignore the error, we will release new bundle soon with workaround.
Customer support service by UserEcho
This is confirmed bug in MATLAB R2019b, see bug-report: https://www.mathworks.com/support/bugreports/details/2088279
I have just updated toolbox bundle for GNU Linux to cope with this mess from TMW. Please re-download toolbox one more time.
Now it includes only proper version of cross, compatible with pre-R2019b and post-R2019b Update 2.
And it is removed from mp.Test to avoid confusion with R2019b.