vpa(4503599627370491.5)
Show older comments
vpa(4503599627370491.5) produces 4503599627370496.0 in 2017b. Why?
Further, sym(4503599627370491.5)-sym(4503599627370496) produces 0.
Answers (5)
Walter Roberson
on 12 Feb 2018
3 votes
sym and vpa are normal functions. Their arguments are evaluated with the normal MATLAB rules before the function is called. Double precision losses information on a number that large before it gets passed on.
You should enclose large numbers in quotes for use with vpa or syms.
Note: the use of sym on character vectors representing numeric constants will continue to be supported, as will sym on a character string representing a simple variable name. The functionality being discontinued within a small number of releases is using sym() on a character string representing an expression. sym() is defined as treating the character strings as MuPad expressions and the replacement functionality str2sym will instead treat it as MATLAB syntax.
Jon Dattorro
on 7 Sep 2023
1 vote
11 Comments
Walter Roberson
on 7 Sep 2023
The MCT user manual contains fundamental misunderstandings of MATLAB, so I would not trust anything it says about MATLAB without further investigation.
This answer doesn't really address the Question. Nevertheless, here is an example from the link above:
format longg
A = rand(3);
A
digits(25)
diff(A)
vpa(A)
diff(vpa(A))
Jon Dattorro
on 7 Sep 2023
Walter Roberson
on 7 Sep 2023
Which continues on to claim that the MATLAB answer is wrong for showing exact zeros. However, vpa(A) is returning sym() objects, and the diff() method of sym() objects is derivative -- so MATLAB is doing exactly what is documented .
The part of the page claiming that diff(vpa(A)) should be returning extended precision numeric differences is displaying a fundamental misunderstanding of how object oriented programming works.
I wrote the company a message of complaint about that section.
Jon Dattorro
on 7 Sep 2023
Pavel Holoborodko
on 19 Sep 2023
Edited: Pavel Holoborodko
on 20 Sep 2023
Dear Walter,
Thank you for your email/complaint. We already replied to your message in private.
But, as the discussion became public, I think it is Ok to publish our reply here.
I cannot publish your email without your permission, so I include only my reply (text between """..."""):
"""
Thank you very much for your comment.
Please be sure, we have a solid understanding of OOP, especially in relation to the MATLAB language. Our toolbox uses OOP extensively, including function overloading, etc. etc.
You are absolutely right about the different meaning of "diff" regarding symbolic expressions. I agree 100%.
However, I cannot agree with you on "expectation" and on how "vpa/diff" should work on numerical inputs.
(1) The original "diff" was one of the first functions introduced in MATLAB with very clear semantics - numeric approximation of the derivatives.
(2) The "VPA" is literally claimed to be "Variable-precision arithmetic (arbitrary-precision arithmetic)" in MATLAB's documentation.
(3) Arbitrary precision is particularly useful for estimating derivatives numerically (as this is a very ill-conditioned problem).
(4) Most vpa-enabled functions check the input type and apply numeric algorithms for numerical inputs. Symbolic inputs treated differently with symbolic-aware algorithms.
Combining the (1)-(4), it is natural to expect that overloaded "vpa/diff" would follow the original meaning and produce a numerical derivative for numerical inputs using extended precision.
The 10+ years of interaction with toolbox users only proves this.
Even prominent researchers, contributing to MATLAB code and close to the team - express this expectation.
In my opinion, to fix this ambiguous situation, the MATLAB team needs to introduce another function for symbolic derivatives (e.g. symdiff).
The "diff" needs to be reserved for numerics, as this semantic has been engraved into users' perception for decades.
***
We will add changes to our documentation regarding the "diff" function (already done).
Thank you,
Pavel.
"""
I can only add, user-accesible documentation for "vpa/diff" method is non-existent. Neither "help vpa.diff", nor "doc vpa.diff" returns meaningful results. The only option for user is to find the actual overloaded code and read the comments. Compare this with function from our toolbox:
>> help mp.diff
diff Difference and approximate derivative.
Usage is identical to built-in function diff.
See also gradient, sum, prod.
To sum up, the sematic of "vpa/diff" function is poorly defined for numeric inputs. This leads to confusion. It has nothing to do with how OOP aspects implemented in MATLAB (we think it is far from pefect, and even published post about it many years ago).
To clear all the confusion, the behavior of "vpa/diff" on numeric inputs must be clearly outlined in the documentation.
Jon Dattorro
on 25 Sep 2023
Please indicate anywhere in any Mathworks documentation that creates a separate class for symbolic numbers that is distinct from sym . In particular, as Pavel expects help vpa/diff to work, please cite anywhere in the documentation that indicates that a class or package named vpa is created.
A = vpa(4503599627370491.5)
class(A)
methods(A)
This shows that vpa() of a numeric value returns a sym, and that it is sym that has method diff
But perhaps it is different for symbolic formulae?
syms x
f = sin(x^2) + gamma(x+1) + 4503599627370491.5
B = vpa(f)
class(B)
C = diff(B)
class(C)
D = vpa(C)
class(D)
Nope, it is sym all the way down, at least in these examples. Is there any other possibility resulting from vpa() ? Only one: in some cases vpa() can return a symfun instead of a sym.
Expecting vpa() to return anything other than sym (or sometimes, symfun) is a substantial misunderstanding of the Symbolic Toolbox
I'm a long time Matlab user and a relatively recent user of the Symbolic Math Toolbox (SMT). Since I started using the SMT, I've had it engraved in my head that SMT diff performs symbolic differentiation whereas base Matlab diff computes numeric differences. I've never confused those two operations. However, having said that, I do recall that when I started using the SMT I did think it odd that TMW chose to use diff for the SMT instead of something that sounds more intuitive, like 'derivative' or 'deriv' or something like that.
What is the function "vpa/diff" that you're trying to find? As I understand it, VPA variables are still of class sym and result in a call to the diff method of the sym class.
x = vpa(pi)
class(x)
which diff(x)
syms x
f(x) = x/sym(pi)
g(x) = vpa(f(x))
Are there any other cases where vpa returns a symfun?
vpa() will return a symfun if the parameter to vpa() evaluates to a symfun .
syms x
f(x) = x/sym(pi)
g1 = vpa(f)
g2 = vpa(3*x + tan(x)^2 + f)
g3 = vpa(f^3)
g4 = vpa( symfun(gamma(x)^2, x) )
g5 = vpa(f(x)) %f(x) is not a symfun
whos
vpa() will not return a symfun in any case where the expression does not evaluate to a symfun. In the example of f(x): invoking a symfun() with a definite parameter evaluates the function on the given parameter, returning a symbolic expression.
Pavel Holoborodko
on 25 Sep 2023
Edited: Pavel Holoborodko
on 25 Sep 2023
@"Expecting vpa() to return anything other than sym (or sometimes, symfun) is a substantial misunderstanding of the Symbolic Toolbox...."
I see, the "vpa" is actually a factory function, generating the instances of "sym" class. I stand corrected, "sym.diff" is the full name of the function.
However, these technical details of SMT implementation have nothing to do with the point we are discussing. We discuss the semantics of the functions in SMT. In other words, what the SMT functions do (or should do) with the input arguments. The type they return is just a technicality in this context.
The SMT does much more than pure symbolic operations. Many functions in SMT apply numerical algorithms (using extended precision) to compute the output results. For example, take a look on EIG:
% LAMBDA = EIG(VPA(A)) and [V,D] = EIG(VPA(A)) compute numeric eigenvalues
% and eigenvectors using variable precision arithmetic. If A does not
% have a full set of eigenvectors, the columns of V will not be linearly
% independent.
Same for SVD, or other cases where result cannot be computed analytically or arbitrary-precision plays an important role.
(Of course the functions still return result of type "sym", but the result itself was computed numerically. I hope this distinction is clear.)
This is exactly the case with "diff". Extended precision is important for computing numerical derivatives/differences accurately, since it is a highly ill-conditioned problem (e.g. same as eigenvalues). Therefore it is quite natural to expect the "sym.diff" to have this functionality.
Besides, the original meaning of "diff" stands the same since the MATLAB inception - computing differences and numeric derivatives. It is not clear why SMT developers decided to change its semantics in SMT. They could've just introduced the new function for symbolic derivatives, as @Paul said. Ambiguity would be avoided.
Here is simple example for illustration:
>> A = vpa(rand(3)); % Store numeric matrix as 'sym'
>> ev = eig(A) % Computes e-vals NUMERICALLY using extended precision (see sym.eig)
ev =
1.752670342542888271728211766393
0.8398632587280311835223044391663
-0.18794383321803341008176836388774
>> class(ev) % Result stored as 'sym', but it is inherently numeric
ans =
'sym'
Similarly, I expect the "sym.diff" to provide the analogous functionality of built-in "diff" and to work as "sym.eig" above - compute result numerically using extended precision. But instead we see:
>> d = diff(A)
d =
[0, 0, 0]
[0, 0, 0]
[0, 0, 0]
11 Comments
James Tursa
on 25 Sep 2023
Edited: James Tursa
on 25 Sep 2023
I don't want this thread to develop into some type of flame war, since I think there are valid points to be made on both sides of the discussion. MATLAB has documented what their sym diff method does, so it is not an error for it to behave as it does. That being said, I am in agreement with @Paul that this was a poor choice on the part of MATLAB and a point of confusion for some that didn't need to be there. Picking "deriv" or something else would have been a better choice. I love working with MATLAB, but I put this one in the bucket of things I wish they had done differently (like making the half class an opaque class instead of a standard numeric class, continuing to hamstring the C-mex interface particularly for classdef objects, etc.).
Walter Roberson
on 25 Sep 2023
Historic note:
MATLAB symbolic diff inherited its name from Maple, from back in the days when the Symbolic Toolbox was mex functions that called into Maplesoft's maple (though I have the suspicion that it might have been from the days of "Waterloo Maple" before they renamed to "Maplesoft")
Still, back then Mathworks could have used a different function name and have rewritten it to diff() on the way to and back from the maple symbolic engine. Might have been a bit tricky to get right in the days when symbolic expressions were mostly quoted strings.
Pavel Holoborodko
on 28 Sep 2023
Apparently the development of extended-precision computations in SMT stuck in the historical times as well:
>> rng(0);
>> digits(25); % Equivalent to 34-digits, as SMT uses 9 guard digits
>> A = vpa(randn(500));
>> tic; [V,D] = eig(A); toc;
Elapsed time is 7016.707970 seconds.
>> tic; [U,S,V] = svd(A); toc;
Elapsed time is 3570.744016 seconds.
These timings are the worst among Maple, Mathematica and MATLAB.
For comparison, here is timings of the (third-party) toolbox mentioned above (MCT):
>> rng(0);
>> mp.Digits(34);
>> A = randn(500,'mp');
>> tic; [V,D] = eig(A); toc;
Elapsed time is 9.198740 seconds.
>> tic; [U,S,V] = svd(A); toc;
Elapsed time is 2.548234 seconds.
* I am the author and main developer of the MCT. These timings can be easily verified since we provide trial versions for tests.
** We have no intentions to use this forum for advertisement or anything alike. The toolbox was mentioned on this thread by another person and it seems legitimate and relevant for discussion to provide more information. Hope this is not against forum rules.
Jon Dattorro
on 28 Sep 2023
Walter Roberson
on 28 Sep 2023
Let me see if I have this straight:
MCT uses a library that slows down MATLAB's double precision arithmetic.
... Therefore MCT should be embedded directly into MATLAB ???
Jon Dattorro
on 29 Sep 2023
Pavel Holoborodko
on 29 Sep 2023
Edited: Pavel Holoborodko
on 29 Sep 2023
Of course, MCT doesn't reduce the speed of double precision computations in MATLAB!
Toolbox overloads standard MATLAB routines for creating elementary arrays - zeros, ones, eye, etc. This is needed to support the 'mp' type as an argument in such functions, e.g. eye(100,'mp'). If requested type is not 'mp' then toolbox dispatches the call to built-in MATLAB functions:
if strcmpi(classname,'mp')
% Create multiprecision array of 'mp' objects
% ...
else
% Call built-in function if standard type is requested
[varargout{1:nargout}] = builtin(funcname,varargin{:}); % <-- 'builtin' is slow
end
As it turns out, the MATLAB's function 'builtin' is really slow. Wrapped call through 'builtin' is ~1.5-2 times slower compared to direct call:
>> tic; for k=1:100000, zeros(100); end; toc; % Direct call
Elapsed time is 0.307563 seconds.
>> tic; for k=1:100000, builtin('zeros',100); end; toc; % Wrapped call
Elapsed time is 0.480524 seconds.
This is where slowdown comes from. It affects only the creation of the elementary arrays, speed of computations are not affected in any way.
Overall, the integration mechanisms provided by MATLAB for third-party toolboxes are quite inefficient and slow. We all know the infamous overhead & limitations of MEX API, especially for custom types (classdef). Slowness of 'builtin' (example above) means that efficient third-party overloads for elementary arrays are (almost) impossible. The list is long.
What @Jon Dattorro is trying to say is that in case of deep integration with MATLAB core, all these inefficiencies can be avoided (for MCT or any other third-party toolbox).
@"Therefore MCT should be embedded directly into MATLAB ???"
Over the last 10+ years MCT has become a de facto standard in the field of extended precision computations in the MATLAB universe, thanks to its quality, speed, breadth and depth of functionality. This could be one of the reasons to consider integrating it into MATLAB's standard set of toolboxes (and deep integration with MATLAB's core).
The core MATLAB functionality is based on the works of third-party developers, starting from Intel MKL, FFTW, UMFPACK/CHOLMOD, MuPAD/Maple and to the ongoing contributions from the researchers from all over the world (e.g. all matrix functions are contributed by Nick Higham team, etc.). The independent researchers/developers are good for the ecosystem, as they are probably much more knowledgeable in the fields they work in.
Steven Lord
on 29 Sep 2023
Toolbox overloads standard MATLAB routines for creating elementary arrays - zeros, ones, eye, etc. This is needed to support the 'mp' type as an argument in such functions, e.g. eye(100,'mp').
If by "overloads" you mean your toolbox defines functions named zeros, ones, eye, etc. that shadow the built-in functions in MATLAB, you should not do that. Follow the instructions given on this documentation page instead to create methods of your class that allow the built-in functions to create instances of your class without affecting the behavior of those built-in functions on other data types.
Pavel Holoborodko
on 29 Sep 2023
Edited: Pavel Holoborodko
on 29 Sep 2023
This is great suggestion, thank you. What version of MATLAB was this introduced? We try to maintain compatibility down to ~R2014-R2016 (our current code with shadowing was written in 2012 and probably compatible downto R2012x).
Is such proper overloading supported for other functions which accept 'classname' argument (e.g. tests matrices - pascal, hadamar, etc.)?
@...without affecting the behavior of those built-in functions on other data types.
We simply call built-in functions for all other data types (except ours). Thus no behavior is affected, besides one extra call to 'builtin'.
Walter Roberson
on 29 Sep 2023
I do not know when it was introduced. I find it documented at least as far back as R2018b, https://www.mathworks.com/help/releases/R2018b/matlab/matlab_oop/class-support-for-array-creation-functions.html (which is as far as the archives go at the moment)
Pavel Holoborodko
on 3 Oct 2023
Edited: Pavel Holoborodko
on 5 Oct 2023
This covers only the simplest case when the user passes the typename argument explicitly (e.g. zeros(...,'mp')).
However, in most of the cases, existing/legacy code relies on default form, without the typename argument.
That is why our users requested to provide a special mode in which the array-creation functions generate extended-precision arrays by default (e.g. when no typename is specified, or it refers to standard floating-point types - 'double' or 'single').
This is useful when legacy code of considerable volume needs to be converted to extended precision:
>> mp.OverrideDoubleBasicArrays(true); % Basic arrays are created as 'mp' by default
>> A = zeros(3);
>> whos A
Name Size Bytes Class Attributes
A 3x3 272 mp
>> % Run legacy code with zeros(...), eye(...), etc.
>> % ...
>> mp.OverrideDoubleBasicArrays(false); % Return to default behavior
That is why our toolbox has to supersede the default functions and call them only when applicable (using 'builtin' calls). This approach works well, with the exception of slow 'builtin' calls.
***
Another interesting case is when the functions are called without arguments at all, e.g. zeros, eye, etc. This case is treated differently by MATLAB - it just returns 0 or 1, without even calling the global/built-in functions. No OOP involved, this is just hardcoded in MATLAB. Meaning these cases are impossible to overload for custom class types at all. Even by superseding the built-in functions.
Jon Dattorro
on 3 Oct 2023
0 votes
Jon Dattorro
on 10 Jan 2025
0 votes
2 Comments
Walter Roberson
on 10 Jan 2025
I have other things on my plate at the moment; I cannot take this on.
Jon Dattorro
on 10 Jan 2025
Categories
Find more on Common Operations in Help Center and File Exchange
Products
on 12 Feb 2018
on 10 Jan 2025
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
