A few days ago, a bug report from our office in Cambridge caught my attention. Computing the singular values and singular vectors of a particular matrix would sometimes cause MATLAB to crash.... read more >>
]]>A few days ago, a bug report from our office in Cambridge caught my attention. Computing the singular values and singular vectors of a particular matrix would sometimes cause MATLAB to crash.
I use two different computers regularly. The machine in my home office is a Lenovo ThinkPad® model T14, loaded with two external monitors, several external disc drives, a sound bar and a dedicated internet connection. My traveling machine is a Lenovo ThinkPad X1 Nano with no external hardware.
The report of a crash in the SVD became even more interesting when I found that it happens on my office computer, but not on the portable. A quick check revealed that the CPUs on the two machines come from different manufacturers. The office computer uses an AMD® Ryzen Pro 5 while the traveling machine uses an Intel® Core i7.
The crash occurs several software layers deep in CGESVD, the LAPACK driver for single precision complex SVD. Various chip manufacturers provide math libraries that have been optimized for their CPUs. However, by default, MATLAB uses the reportedly faster Intel Math Kernel Library, MKL. It is possible to switch to other math libraries.
We have experts at MathWorks who know far more about the details of these libraries than I do. They will soon sort this all out. In the meantime, here is what I have learned about the offending matrix.
We refer to the matrix in the crash report by its case number in our bug tracking system. The matrix is of modest size but is otherwise unusual for several reasons. It is rectangular with fewer rows than columns, it is in single precision, and it is complex.
clear load g3366394 X whos
Name Size Bytes Class Attributes X 219x384 672768 single complex
The following code calling SVD with three outputs will crash MATLAB on my T14, but not on my X1.
trysvd = false if trysvd [U,S,V] = svd(X); R = U*S*V' - X; relres = norm(R)/norm(X) end
trysvd = logical 0
Computing the singular values without the vectors can be done on either machine. The following code uses double precision and then plots the singular values on a logarithmic scale with a line at single precision roundoff level.
S = svd(X); semilogy(S,'.-') ep = eps('single'); line([0 230],[ep ep]) axis([0 230 1e-15 10]) legend({'singular values','eps(''single'')'})
We see that the matrix is far from full rank. About half of its singular values are negligible. This is confirmed by
xrank = rank(X)
xrank = 110
The cause of the low rank is easy to find. This surf plot reveals that almost half of the rows are flat zero.
Let's remove the zero rows.
c = sum(abs(X),2)==0; nnzc = nnz(c) X(c>0,:) = [];
nnzc = 109
The remaining matrix is full rank and it is safe to compute its singular vectors.
[U,S,V] = svd(X); R = U*S*V' - X; resnorm = norm(R)
resnorm = 2.8191e-06
Removing the zero rows was the first work-around that I tried. There are many others. You can replace the zeros with any nonzero "fuzz".
fuzz = 1.e-20; [U,S,V] = svd(X + fuzz*randn(size(X)));
resnorm = norm(U*S*V'-X)
resnorm = 2.8222e-06
You can reorder the matrix so that the zero rows are not at the beginning.
[U,S,V] = svd(flipud(X)); U = flipud(U);
resnorm = norm(U*S*V'-X)
resnorm = 2.3809e-06
How to avoid the crash is not the most important question. What causes the crash with the original matrix? We will find out and get it fixed.
Get
the MATLAB code
Published with MATLAB® R2024a
Floating point arithmetic doesn't get the respect it deserves. Many people consider it mysterious, fuzzy, unpredictable. These misgivings often occur in discussion of vector sums. Our provocatively named SuperSum is intended to calm these fears.... read more >>
]]>Floating point arithmetic doesn't get the respect it deserves. Many people consider it mysterious, fuzzy, unpredictable. These misgivings often occur in discussion of vector sums. Our provocatively named SuperSum is intended to calm these fears.
A ledger is a list of transactions in an account. Auditing the ledger involves comparing the total of the items with the change in the account balance.
If v is a vector of transaction amounts, then we need to compute
sum(v)
If this sum is equal to the change in the balance, then it is reasonable to assume that the ledger is correct. If not, the ledger must be examined line-by-line.
Have you ever used a checksum for a file transfer? One checksum is computed before the file is sent. After the file has been sent over a questionable channel, a second checksum is computed on the receiving end. If the two checksums agree, the transmission was probably okay. If not, the file must be sent again.
Floating point addition is not associative. This means
(a + b) + c
is not necessarily the same as
a + (b + c).
So the order of doing a computation is important.
Computers are deterministic devices. If the same computation is done in the same order more than once, the results must be the same. If you get different sums from different runs on a fixed computer, then it must be because the order has been changed.
In recent years, we have made the built-in function sum(x) faster by parallelizing it. The input vector is broken into several pieces, the sum of each piece is computed separately and simultaneously, then the partial sums are combined to give the final result. If the number and size of the pieces varies from run to run, the order varies and consequently the computed sums may vary.
Here are three replacements for nansum(x), the version of sum(x) that skips over NaNs and infs.
simplesum
Avoid the effect of blocking in the built-in sum(x).
function s = simplesum(x) % simplesum. s = simplesum(x), for vector(x). % Force sequential order for sum(x). % Skips NaNs and infs.
s = 0; for k = 1:length(x) if isfinite(x(k)) s = s + x(k); end end end
KahanSum
This is the Kahan summation algorithm. The sum is accumulated in two words, the more significant s and the correction c. If you were able to form s + c exactly, it would be more accurate than s alone.
function s = KahanSum(x) % KahanSum. s = KahanSum(x), for vector(x). % More accurate, but slower, than sum(x). % Skips NaNs and infs. % https://en.wikipedia.org/wiki/Kahan_summation_algorithm.
s = 0; % sum c = 0; % correction for k = 1:length(x) if isfinite(x(k)) y = x(k) - c; t = s + y; c = (t - s) - y; s = t; end end end
SuperSum
I gave it a pretentious name to grab attention. Use extended precision, with enough digits to hold any MATLAB double.
function s = SuperSum(x) % SuperSum. s = SuperSum(x), for vector(x). % Symbolic Math Toolbox extended precision. % Skips NaNs and infs. % % 632 decimal digits holds every IEEE-754 double. % 632 = ceil(log10(realmax) - log10(eps*realmin)); % din = digits(632); % Remember current setting s = double(sum(sym(x(isfinite(x)),'D'))); digits(din) % Restore end
A test case here at MathWorks, known by a French-Canadian name that translates to "toy example", has a vector e2 of length 4243 and values that range from -3.3e7 to 6.8e9.
When running tests involving floating point numbers it is a good idea to set the output format to hex so we can see every last bit.
format hex load jeuTest e2 x = e2;
[nansum(x) simplesum(x) KahanSum(x) SuperSum(x)]
ans = 423785e8206150e2 423785e8206150e0 423785e8206150e1 423785e8206150e1
For jeuTest, Kahan summation gets the same result as SuperSum, while nansum and simplesum differ in the last bit or two.
The vector length is only three, but the third term completely cancels the first, and the second term rises from obscurity. In this situation, KahanSum is no more accurate than sum.
format hex x = [1 1e-14 -1]'
[nansum(x) simplesum(x) KahanSum(x) SuperSum(x)]
x = 3ff0000000000000 3d06849b86a12b9b bff0000000000000
ans = 3d06800000000000 3d06800000000000 3d06800000000000 3d06849b86a12b9b
I will leave careful timing for another day. Let's just say that in situations like jeuTest, KahanSum is probably all you need. It is usually more accurate than sum, and not much slower.
However, for complete reliability, use SuperSum. There is no substitute for the right answer.
Get
the MATLAB code
Published with MATLAB® R2024a
Our technical support group recently received a request for a tool that would convert IBM System/360 hexadecimal floating point numbers to the IEEE-754 format. I am probably the only one left at MathWorks that actually used IBM mainframe computers. I thought we had seen the last of hexadecimal arithmetic years ago. But, it turns out that the hexadecimal floating point format is alive and well.... read more >>
]]>Our technical support group recently received a request for a tool that would convert IBM System/360 hexadecimal floating point numbers to the IEEE-754 format. I am probably the only one left at MathWorks that actually used IBM mainframe computers. I thought we had seen the last of hexadecimal arithmetic years ago. But, it turns out that the hexadecimal floating point format is alive and well.
The System/360 is a family of mainframe computers that IBM introduced in 1965 and that dominated the computer industry until PCs came along twenty years later. They range in size from desk-sized to systems that fill a large room.
Here is a photo of a mid-sized model.
System/360, Model 60. Photo from Ken Shirrif's blog, IBM 360/System Summary.
The System/360 architecture is byte-oriented, so it can handle business data processing as well as scientific and engineering computation. This leads to base-16, rather than base-2 or base-10, floating point arithmetic.
* Binary f*2^e 1/2<=f<1 * Decimal f*10^e 1/10<=f<1 * Hexadecimal f*16^e 1/16<=f<1
Floating point formats played an important role in technical computing in the early days. This table from FMM lists formats that were in use in the 1970s, before IEEE-754 was introduced in 1985.
The System/360 hexadecimal format is used in many industries for the preservation of data files.
CREWES. Teaching exploration seismology. Comprehensive MATLAB toolbox for use with the textbook "Numerical Methods of Exploration Seismology with algorithms in Matlab" by Gary F. Margrave, a geoscience professor at the University of Calgary.
Library of Congress. Government.
NASA. Astronautics.
SAS. Statistics and business analytics. SAS wrapers for C.
Enthought. Python wrappers for C.
Hex_ieee. I have two twenty-line MATLAB functions, ieee2ibm and ibm2ieee, that convert IEEE-754 floating point to and from IBM hexadecimal format.
Three statements in the middle of ieee2ibm are the key to the entire operation. The first statement is
[~,e] = log2(x)
With two output arguments, log2 returns the mantissa and exponent of an IEEE-754 floating point number. The mantissa is not needed here.
The second key statement
e = ceil(e/4)
makes e divisible by 4. This turns e into the appropriate hexadecimal exponent so that the third statement
f = x.*16.^(-e)
can produce the hexadecimal mantissa.
function z = ieee2ibm(x) Convert IEEE-754 to IBM System 360 hexadecimal. z = ieee2ibm(x) Input x, real column vector. Output z, length(x)-by-16 char. Example: ieee2ibm(-118.625) = 'C276A00000000000'.
s = sign(x); % -1, 0, or 1 x = abs(x); x(x < 16^(-65)) = 0; % Underflow x(x >= 16^63) = (1-eps/2)*16^63; % Overflow
[~,e] = log2(x); % base 2 exponent e = ceil(e/4) % base 16 exponent f = x.*16.^(-e); % base 16 mantissa
E = uint64((e+64)*2^56); % Assemb1e output F = uint64(f*2^56); S = uint64((1-s)*2^62); % 1 or 0 z = dec2hex(S + E + F); % z = 'ZZFFFFFFFFFFFFFF' end
function x = ibm2ieee(z) Convert IBM System 360 hexadecimal to IEEE-754. x = ibm2ieee(z) Input z, n-by-16 char. Output x, n-by-1 double. Example: ibm2ieee('C276A00000000000') = -118.625.
E = hex2dec(z(:,1:2)); % Disassemble input F1 = hex2dec(z(:,3:8)); % < 16^6 F2 = hex2dec(z(:,9:end)); % < 16^8 s = sign(128-E); % -1 or 1
e = E-(s>0)*64-(s<0)*192; % base 16 exponent f = F1/16^6 + F2/16^14; % base 16 mantissa x = s.*f.*16.^e; end
Underflow. Anything < 16^(-65) is too small and is flushed to zero. There are no denormals.
Overflow. Anything >= 16^63 is too large. There is no inf or NaN.
* 1.0 4110000000000000 * 0.1 401999999999999A * -pi C13243F6A8885A30 * 5.3976e-79 0010000000000000 * 7.2370e+75 7FFFFFFFFFFFFFF8
S/360 hexadecimal has 7 exponent bits, while IEEE-754 has 11. Consequently, hexadecimal has a much smaller range, 5.4e-79 to 7.2e+75 versus 2.2e-308 to 1.8e+308.
The base-16 normalization implies that hexadecimal effectively has between 53 and 56 mantissa bits. Counting the hidden bit, IEEE-754 also has 53. So, the accuracy of the two is pretty much the same.
My functions ieee2ibm and ieee2ibm described above, modified to handle both single and double, plus hex_test, which does what its name implies, are available at Hex_ieee.
Homework: What happens?
ok = 0; for k = 1:10 x = single(k/10); ok(k) = hex_test(x); end ok
Get
the MATLAB code
Published with MATLAB® R2024a
The graphics in my post about R^2 were produced by an updated version of a sixty-year old program involving the U.S. census. Originally, the program was based on census data from 1900 to 1960 and sought to predict the population in 1970. The software back then was written in Fortran, the predominate technical programming language a half century ago. I have updated the MATLAB version of the program so that it now uses census data from 1900 to 2020.... read more >>
]]>The graphics in my post about R^2 were produced by an updated version of a sixty-year old program involving the U.S. census. Originally, the program was based on census data from 1900 to 1960 and sought to predict the population in 1970. The software back then was written in Fortran, the predominate technical programming language a half century ago. I have updated the MATLAB version of the program so that it now uses census data from 1900 to 2020.
The latest version of the census application is now available at censusapp2024. Here are the data and the opening screenshot.
[t,p]=UScensus;fprintf('%12d%12.3f\n',[t,p]')
1900 75.995 1910 91.972 1920 105.711 1930 123.203 1940 131.669 1950 150.697 1960 179.323 1970 203.212 1980 226.505 1990 249.633 2000 281.422 2010 308.746 2020 331.449
Today, MATLAB makes it easier to vary parameters and visualize results, but the underlying mathematical principles are unchanged:
One new observation is added to the data every 10 years, when the United States does the decennial census. Originally there were only 7 observations; today there are 13. The program now allows you to fit the data exactly by interpolation with a polynomial of degree 12 or fit it approximately by polynomials of degree less than 12.
Here are the least-squares fits with linear, cubic, and degree seven polynomials and the interpolating polynomial. As the polynomial degree increases, so does R^2, until R^2 reaches one with the exact fit.
Do any of these fits look like they could be used to predict future population growth?
In addition to polynomials, you can choose interpolation by three different piecewise Hermite cubics.
Since these fits interpolate the data, all their R^2 values are one. But before 1900 and after 2020 these functions are cubic polynomials that are not designed for extrapolation.
It is also possible to do nonlinear least squares fits by an exponential, a logistic sigmoid, and an exponential of an exponetial known as the Gompertz model.
An article by Kathleen and Even Tjørve, from the Inland Norway University of Applied Sciences in Elverum, Norway, in the journal PLOS ONE has this to say about Gompertz. "The Gompertz model has been in use as a growth model even longer than its better known relative, the logistic model. The model, referred to at the time as the Gompertz theoretical law of mortality, was first suggested and first applied by Mr. Benjamin Gompertz in 1825. He fitted it to the relationship between increasing death rate and age, what he referred to as 'the average exhaustions of a man’s power to avoid death” or the 'portion of his remaining power to oppose destruction.' "
Which fits are suitable for predicting future population size?
Despite their large R^2 values, polynomials of any degree are not suitable because outside of the time interval they behave like polynomials and do not provide realistic predictions.
Splines were never intended for extrapolation.
That leaves the exponentials. The simple exponential model grows exponentially and is not suitable. The Gompertz fit does approach a finite asymptotic limit, but the value is an astronimical a = 2101, corresponding to 2.1 $\times 10^9$ inhabitants. Hopefully, that is out of the question.
The logistic fit has an asymptotic limit of a = 655.7. We recently passed the value of t where p(t) reaches a/2, namely c = 2018. So, the logistic model predicts that the long-term size of the U.S. population will be about twice its current value. Is that realistic? Probably not.
The British statistician George Box once said, "all models are wrong, some are useful." This is true of the models of the U. S. Census that I have discussed over the past sixty years.
Here is censusapp2024 after all its buttons have been pushed. The extrapolation date is set to 2040. White noise has been added to the data. The model is a fourth-degree polynomial with an R^2 = 0.99. The R^2 value and the error estimates produced by errs account for errors in the data, but not in the model.
This particular model does a lousy job of predicting even twenty years in the future. Some of the other models are better, many are worse. Hopefully, their study is worthwhile.
I have made blog posts about the census before, in 2020 and in 2017.
Predicting population growth is featured in Computer Methods for Mathematical Computations, by George Forsythe, Mike Malcolm and myself, published by Prentice-Hall in 1977. That textbook is now available from an interesting smorgasbord of sources, including Google Scholar, Amazon, dizhasneatstuff, Abe Books, Internet Archive, PDAS, WorldCat (Chinese).
censusapp2024 is available at censusapp2024.
Get
the MATLAB code
Published with MATLAB® R2024a
The coefficient of determination, R-squared or R^2, is a popular statistic that describes how well a regression model fits data. It measures the proportion of variation in data that is predicted by a model. However, that is all that R^2 measures. It is not appropriate for any other use. For example, it does not support extrapolation beyond the domain of the data. It does not suggest that one model is preferable to another.... read more >>
]]>The coefficient of determination, R-squared or R^2, is a popular statistic that describes how well a regression model fits data. It measures the proportion of variation in data that is predicted by a model. However, that is all that R^2 measures. It is not appropriate for any other use. For example, it does not support extrapolation beyond the domain of the data. It does not suggest that one model is preferable to another.
I recently watched high school students participate in the final round of a national mathematical modeling competition. The teams' presentations were excellent; they were well-prepared, mathematically sophisticated, and informative. Unfortunately, many of the presentations abused R^2. It was used to compare different fits, to justify extrapolation, and to recommend public policy.
This was not the first time that I have seen abuses of R^2. As educators and authors of mathematical software, we must do more to expose its limitations. There are dozens of pages and videos on the web describing R^2, but few of them warn about possible misuse.
R^2 is easily computed. If y is a vector of observations, f is a fit to the data and ybar = mean(y), then
R^2 = 1 - norm(y-f)^2/norm(y-ybar)^2
If the data are centered, then ybar = 0 and R^2 is between zero and one.
One of my favorite examples is the United States Census. Here is the population, in millions, every ten years since 1900.
t p ____ _______ 1900 75.995 1910 91.972 1920 105.711 1930 123.203 1940 131.669 1950 150.697 1960 179.323 1970 203.212 1980 226.505 1990 249.633 2000 281.422 2010 308.746 2020 331.449
There are 13 observations. So, we can do a least-squares fit by a polynomial of any degree less than 12 and can interpolate by a polynomial of degree 12. Here are four such fits and the corresponding R^2 values. As the degree increases, so does R^2. Interpolation fits the data exactly and earns a perfect core.
Which fit would you choose to predict the population in 2030, or even to estimate the population between census years?
R2_census
Thanks to Peter Perkins and Tom Lane for help with this post.
Get
the MATLAB code
Published with MATLAB® R2024a
The Closest Pair of Points problem is a standard topic in an algorithms course today, but when I taught such a course fifty years ago, the algorithm was not yet known.... read more >>
]]>The Closest Pair of Points problem is a standard topic in an algorithms course today, but when I taught such a course fifty years ago, the algorithm was not yet known.
Imagine you are driving a car on the Harbor Freeway in southern California with typical Los Angeles traffic conditions. Among the many things you might want to know is which pair of vehicles is nearest each other.
This is an instance of the Closest Pair of Points problem:
It is convenient to represent the points by a vector of complex values. The distance between points z(k) and z(j) is then
d = abs(z(k) - z(j))
Here are a few points in the unit square. The closest pair is highlighted.
The first algorithm you might think of computes the distances between all possible pairs of points and finds the minimum. This is a brute force approach that requires only a few lines of code.
function d = Pairs(z) % Pairs. % d = Pairs(z) is the minimum distance between any two elements % of the complex vector z.
n = length(z); d = Inf; for k = 1:n for j = k+1:n if abs(z(k) - z(j)) < d d = abs(z(k) - z(j)); end end end end
DivCon stands for Divide and Conquer. In outline, the steps are:
function d = DivCon(z,sorted) % DivCon. % d = DivCon(z) is the minimum distance between any two elements % of the complex vector z. % % d = DivCon(z,true) is a recursive call with ascending real(z).
n = length(z); if n <= 3 d = Pairs(z); else if nargin < 2 || ~sorted [~,p] = sort(real(z)); z = z(p); end m = floor(n/2);
% Left half dl = DivCon(z(1:m),true)
% Right half dr = DivCon(z(m+1:end),true);
% Choose d = min(dl,dr);
% Center strip ds = Center(z,d); d = min(ds,d); end end
The delicate case involves the strip of points near the center dividing line. The width of the strip is the closest distance found in the recursion. Any closer pair with one point in each half must be in this strip.
function d = Center(z,d) % Center(z,d) is used by DivCon to examine the % strip of half-width d about the center point.
n = length(z) m = floor(n/2); xh = real(z(m)); [~,p] = sort(imag(z)); z = z(p); s = []; for i = 1:n if abs(real(z(i)) - xh) < d s = [s; z(i)]; end end
ns = length(s); for k = 1:ns for j = k+1:ns if (imag(s(j)) - imag(s(k))) < d && abs(s(k) - s(j)) < d d = abs(s(k) - s(j)); end end end end
Let n be the number of points. An asymptotic execution-time complexity analysis involves n approaching infinity.
It is not hard to see that the complexity of the brute force algorithm implemented in Pairs is O(n^2).
There are dozens of pages on the web devoted to showing that the complexity of the divide and conquer algorithm implemented in DivCon and Center is O(n*log(n)). The best page that I have seen is the YouTube video by Ling Qi. The key to the analysis is showing that the inner loop in Center is executed at most 7 times for any n.
We measured the execution time of Pairs(z) and DivCon(z) for n from 1,000 to 40,000 and computed the ratios of the two times. The complexity analysis predicts that this ratio is asymptotically
O(n/log(n))
Here are the timing results and a least square fit by n/log(n).
A self-extracting MATLAB archive is available at https://blogs.mathworks.com/cleve/files/TestDivCon_mzip.m
Ling Qi, IDeer7, Closest Pair of Points (Divide and Conquer) Explained. https://www.youtube.com/watch?v=6u_hWxbOc7E.
Cormen, Thomas H.; Leiserson, Charles E.; Rivest, Ronald L.; Stein, Clifford. Introduction to Algorithms (4th ed.). MIT Press and McGraw-Hill. ISBN 0-262-04630-X. 1312 pp.
Get
the MATLAB code
Published with MATLAB® R2023a
I have just returned from the MathWorks company meeting celebrating our 40th Anniversary. In one of the presentations, Jos Martin described how Parallel MATLAB was introduced almost twenty years ago. Here are a few slides from Jos's talk.... read more >>
]]>I have just returned from the MathWorks company meeting celebrating our 40th Anniversary. In one of the presentations, Jos Martin described how Parallel MATLAB was introduced almost twenty years ago. Here are a few slides from Jos's talk.
In MATLAB News and Notes for spring 1995, I wrote a one-page Cleve's Corner titled "Why there isn't any parallel MATLAB." There were three reasons.
This one-page note turned out to be one of my most widely cited publications.
.
.
A 2001 survey by Ron Choy at MIT found 27 different projects that were developing some way to run MATLAB in parallel. All of them involved a MATLAB-based host program calling a fixed library of parallel functions, written in some other language, on the workers. None of the systems were capable of running arbitrary MATLAB programs in parallel. None of them were MathWorks products.
MathWorks introduced the MATLAB Distributed Computing Toolbox in November 2004. We improvised this demo setup at our first Supercomputing Conference, SC2004 in Pittsburg,
.
A year later, SC2005 was in Seattle and our booth featured four worker machines on a wire shelving unit purchased at a local Home Depot.
.
Since Seattle was his home town, Bill Gates gave the keynote talk at SC2005. He announced that Microsoft was going into High Performance Computing and used the MathWorks Distributed Computing Toolbox in his demonstration.
.
So, a little more than ten years after the first Cleve's Corner about parallel computing, a second Cleve's Corner in News and Notes was able to reverse the situation.
.
Get
the MATLAB code
Published with MATLAB® R2023a
The Swinging Sticks is a kinetic sculpture that exhibits chaotic motion. The device became very popular after it upstaged Tony Stark in Iron Man 2. My daughter Carolyn gave me a desktop version of Swinging Sticks for Christmas. I immediately set out to simulate it.... read more >>
]]>The Swinging Sticks is a kinetic sculpture that exhibits chaotic motion. The device became very popular after it upstaged Tony Stark in Iron Man 2. My daughter Carolyn gave me a desktop version of Swinging Sticks for Christmas. I immediately set out to simulate it.
Chaotic motion appears random but isn't. Once the motion begins, the initial conditions together with Newton's law of motion, F = ma, determine subsequent behavior. There are no random forces. It may be difficult to predict positions, but they are well-determined, nonetheless.
A classic example of chaotic motion is the double pendulum. One mass at end of a massless string swings about a fixed pivot, and a second mass is attached by a massless string to the first. My simulator of the classic double pendulum is available in Cleve's Lab, swinger, and a movie is available here pendulum movie.
The swinging sticks are similar to the double pendulum. The sticks are two rods with uniformly distributed mass, different lengths and off-center pivots. The best way to view the motion is to download this code and run it in your own MATLAB. Otherwise, here in a short slow-motion animated GIF.
And, here is a longer Swinging Sticks Video.
The motion of the shorter of the two rods is chaotic. Here are the orbits traced by the ends of the short rod.
Swinging Sticks sculptures are available in various sizes and colors. The Swinging Sticks.
Our mathematical model is of a frictionless perpetual motion machine. The real sculptures have an ingenious electromagnetic controller in the base that is claimed to run for two years on four AA batteries. Mine has been running since Christmas. An excellent YouTube video by Wayne Schmidt describes the controller.
https://blogs.mathworks.com/cleve/files/swinging_sticks.m
Get
the MATLAB code
Published with MATLAB® R2024a
Nick Higham passed away last Saturday. Nick was a close friend of mine and a great friend of MATLAB. I will leave it to others to describe his research and teaching, his many honors, and his service to our community, especially SIAM. I have just a few, more personal, comments.... read more >>
]]>Nick Higham passed away last Saturday. Nick was a close friend of mine and a great friend of MATLAB. I will leave it to others to describe his research and teaching, his many honors, and his service to our community, especially SIAM. I have just a few, more personal, comments.
Monday's NA Digest led off with this from Nick's wife Francoise and his brother Des.
Subject: Nick Higham (1961--2024)
With great sadness we report that Nick Higham, Royal Society Research Professor and Richardson Professor of Applied Mathematics at the University of Manchester, passed away on January 20, 2024, at the age of 62 after an 18 month struggle with a form of blood cancer. An obituary describing Nick's research and leadership contributions will appear in SIAM News in due course.
Francoise Tisseur and Des Higham
Nick was an excellent writer, and an excellent writer about writing.
Here are the covers of his six books.
SIAM published five of these. Two are surveys of Nick's research on the accuracy of numeric algorithms and the computation of matrix functions. Two more, one of them coauthored with Dennis Sherwood, are guides to mathematical exposition.
MATLAB Guide, by Des and Nick Higham, is one of my favorite books about MATLAB. It is a succinct introduction for newcomers and a valuable refresher for old-timers. The third edition, published in 2017, includes chapters on object-oriented computing, parallel computing, the Symbolic Math Toolbox and other recent additions. Be sure to check out the MATLAB Guide web site.
The only non-SIAM book pictured above is The Princeton Companion to Applied Mathematics. It is over 1,000 pages long and features nearly 200 sections written by an international team of experts. Nick is the editor-in-chief and wrote many of the sections himself.
Here is a Word Cloud from Nick's home page. It shows the frequency of the tags for his blog and confirms his interest in MATLAB.
Anyone interested in numerical linear algebra should also be interested in the gallery function, which is based on Nick's work. Enter
>> doc gallery
Scroll down to matrixname and investigate over 70 different test matrices.
If you find gallery irresistible, take a look at anymatrix, an extensible matrix collection, by Nick and Mantas Mikaitis.
This is very personal for me. Thirty or forty years ago, Charlie Van Loan and I were regarded as authorities on computing the matrix exponential, $e^{A}$. The function expm has been in MATLAB since its very beginning. Around twenty years ago, we ceded the authority title to Nick and Awad Al-Mohy. Their code for matrix exponential is now the basis for expm.
Our business has lost one of its superstars. I have lost a good friend, way too soon. Goodbye Nick.
Get
the MATLAB code
Published with MATLAB® R2023a
We have been investigating a recent bug report about fitnlm, the Statistics and Machine Learning Toolbox function for robust fitting of nonlinear models.... read more >>
]]>We have been investigating a recent bug report about fitnlm, the Statistics and Machine Learning Toolbox function for robust fitting of nonlinear models.
The bug report comes from Greg Pelletier, an independent research scientist and biogeochemical modeler in Olympia, Washington. Greg has been studying the vulnerability of sensitive marine organisms to increases in ocean acidification. One of the most important of these organisms is Mercenaria mercenaria, the hard clam.
Especially here in New England, hard clams are known by their traditional Native American name, quahog. They have a well-deserved reputation for making excellent clam chowder.
We are all aware of increasing levels of carbon dioxide in the earth's atmosphere. We may not be as aware of the effect this increase has on the health of the earth's oceans. According to NOAA, the ocean absorbs about 30% of the atmospheric carbon dioxide.
A definitive and controversial 2009 paper by Justin Ries and colleagues, then at the Woods Hole Oceanographic Institution, is "Marine calcifiers exhibit mixed responses to CO2-induced ocean acidification", https://doi.org/10.1130/G30210A.1. The hard clam example in Greg's bug report comes from figure 1K in the Ries et al. paper.
The independent variable in experiments is the ratio of alkalinity of sea water to the concentration of dissolved inorganic carbon. The dependent variable is the calcification rate, which compares how fast the organism builds its shells to how fast the shells are dissolving.
The model chosen by Ries at al. is
$$ y \approx \beta_1 + \beta_2 e^{\lambda t} $$
where $t$ is the ratio of alkalinity to dissolved carbon and $y$ is the calcification rate. The data have only four distinct values of $t$, with several observations of $y$ at each value.
The parameters $\beta_1$, $\beta_2$ and $\lambda$ are determined by least squares curve fit. This is a separable least squares problem. For any given value of $\lambda$, the parameters $\beta_1$ and $\beta_2$ occur linearly and the least squares solution can be obtained by MATLAB's backslash.
Gene Golub and Victor Pereyra described separable least squares in 1973 and proposed solving it by a variable projection algorithm. Since 1973 a number of people, including Pereyra, Linda Kaufman, Fred Krogh, John Bolstadt and David Gay, have contributed to the development of a series of Fortran programs named varpro. In 2011, Dianne O'Leary and Burt Rust created a MATLAB version of varpro. Their report, https://www.cs.umd.edu/~oleary/software/varpro/, is a good background source, as well as documentation for varpro.m.
I have a section on separable least squares, and an example, expfitdemo, in NCM, Numerical Computing with MATLAB. I have modified expfitdemo to work on Greg's quahogs problem.
It turns out that the problem Greg encountered can be traced to the fact that the data are not centered. The given values of $t$ are all positive. This causes fitnlm to print a warning message and attempt to rectify the situation by changing the degrees of freedom from 22 to 23, but this only makes the situation worse. (We should take another look at the portion of fitnlm that adjusts the degrees of freedom.)
It is always a good idea in curve fitting to center the data with something like
t = t - mean(t)
The values of $y$ are already pretty well centered. Rescaling $y$ with
y = 10000*y
makes interpretation of results easier.
With the data centered and scaled, we have three different ways of tackling Greg's problem. All three methods agree on the results they compute.
Nonlinear regression model: y ~ param1 + param2*exp(param3*xval)
Estimated Coefficients: Estimate SE tStat pValue ________ _______ _______ __________
param1 0.69536 0.1657 4.1964 0.00037344 param2 -0.26482 0.19909 -1.3302 0.19709 param3 -22.218 8.1494 -2.7263 0.012327
Number of observations: 25, Error degrees of freedom: 22 Root Mean Squared Error: 0.307 R-Squared: 0.828, Adjusted R-Squared 0.813 F-statistic vs. constant model: 53, p-value = 3.86e-09
Linear Parameters: 0.695367 -0.264837 Nonlinear Parameters: -22.217495
Norm of weighted residual = 1.438935 Norm of data vector = 3.545820 Expected error of observations = 0.306782 Coefficient of determination = 0.828145
Regression.t_ratio 4.1962 -1.3301 -2.7264
Regression.std_param 0.1657 0.1991 8.1490
lambda = -22.2180 condX = 4.3882 beta = 0.6954 -0.2648 normres = 1.4389
quahogsfit produces this plot, which can be compared with figure 1K from Ries et al, reproduced above.
The codes for this post are available here quahogs_driver.m and here varpro.m.
Thanks to Greg Pelletier for the bug report and to Tom Lane for his statistical expertise.
Get
the MATLAB code
Published with MATLAB® R2023a