Thank you. The coefficients c returned by J are too off so I did not check the absolute error.
Yes, the absolute error is given by >./|(X (+/ .*) c)-^x The result is: x=:_0.5+(i.1000)%999 c=:(^x) %. X=:(^ 0j1 * x */ i:10) >./|(X (+/ .*) c)-^x 0.683384 The least squares method minimizes not the absolute error but the mean squared error, or, (+/%#)*:|(X (+/ .*) c)-(^x) 0.0290534 The corresponding error in Octave is >> mean(abs(X*c-exp(x)).^2) ans = 1.8285e-30 >> Thanks. Yours sincerely, Imre Patyi On Thu, Apr 29, 2021 at 4:17 PM Raul Miller <[email protected]> wrote: > I do not see a test expression in the J implementation. > > Is it fair to say that the octave expression > > norm(X*c-exp(x),"inf") > > represents the value which should be minimized? > > If so, does this J expression match what you were thinking of? > +/|(X+/ .*c)-^x > > Thanks, > > -- > Raul > > On Thu, Apr 29, 2021 at 12:17 AM Imre Patyi <[email protected]> wrote: > > > > Dear [email protected] , > > > > I have a few years' worth of toy experience with J. > > I was playing with the toy problem of fitting a trigonometric > > polynomial with complex coefficients to the curve y=e^x, > > -1/2 <=x<= 1/2 through 1000 sample points. > > I tried using the linear least squares command %. but it > > seems to me that it gives coefficients c several orders of > > magnitude larger than what seems right. > > > > Could you please tell me what I would need to change to fix the > > following little J program (e.g., how to replace %. with a call > > to LAPACK)? (I am thinking that J might not have a > > numerically stable QR decomposition or SVD behind %. --- I tried to > > find the source code file for %. and only found > > > > https://github.com/jsoftware/jsource/blob/master/jsrc/vd.c > > > > which talks about qr decomposition and a function qrr(), but > > I did not spend too much time on trying to locate the source code > > for qrr() ). > > > > You can see from the comparable Octave and Python versions that > > the result can be much more accurate than what I got out of J > > for this ill-conditioned problem. I also tried it on a raspberry > > pi with a newer version of J, but it was still several orders > > of magnitude wrong. > > > > Thank you very much! > > > > Yours sincerely, > > Imre Patyi > > > > ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ > > > > x=:_0.5+(i.1000)%999 > > c=:(^x) %. X=:(^ 0j1 * x */ i:10) > > (i.5){"1(i.4){X > > ,.c > > > > NB. Sample run. Windows 10, i5, J805. > > NB. x=:_0.5+(i.1000)%999 > > NB. c=:(^x) %. X=:(^ 0j1 * x */ i:10) > > NB. (i.5){"1(i.4){X > > NB. 0.283662j_0.958924 _0.210796j_0.97753 _0.653644j_0.756802 > > _0.936457j_0.350783 _0.989992j0.14112 > > NB. 0.274049j_0.961716 _0.219594j_0.975591 _0.659683j_0.751544 > > _0.938892j_0.344213 _0.989127j0.147063 > > NB. 0.264409j_0.964411 _0.228374j_0.973574 _0.66568j_0.746237 > > _0.94128j_0.337626 _0.988226j0.153001 > > NB. 0.254742j_0.967009 _0.237135j_0.971477 _0.671635j_0.740882 > > _0.943623j_0.331022 _0.987289j0.158934 > > NB. ,.c > > NB. 4.05363e6j_2.63539e6 > > NB. _5.04244e7j3.27836e7 > > NB. 2.90519e8j_1.88887e8 > > NB. _1.02273e9j6.64963e8 > > NB. 2.44111e9j_1.58715e9 > > NB. _4.13678e9j2.6895e9 > > NB. 5.05136e9j_3.2837e9 > > NB. _4.39569e9j2.8568e9 > > NB. 2.60811e9j_1.69431e9 > > NB. _9.47205e8j6.14882e8 > > NB. 1.5859e8j_1.02848e8 > > NB. _1.84539e7j1.19961e7 > > NB. 7.13468e7j_4.62898e7 > > NB. _1.82918e8j1.19028e8 > > NB. 3.31731e8j_2.16376e8 > > NB. _3.96625e8j2.58975e8 > > NB. 3.10552e8j_2.02872e8 > > NB. _1.60195e8j1.04678e8 > > NB. 5.30711e7j_3.46859e7 > > NB. _1.03281e7j6.75135e6 > > NB. 904018j_591051 > > > > ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ > > > > x=linspace(-0.5,0.5,1000)' ; X=exp(1i*x.*(-10:10)) ; > > X(1:4,1:5) > > c=X\exp(x) > > norm(X*c-exp(x),"inf") > > f=@(x)exp(1i*x.*(-10:10))*c ; > > y=linspace(-0.5,0.5,1e6)';norm(f(y)-exp(y),"inf") > > > > %{ > > Sample run. > > GNU Octave, version 5.2.0 > > Copyright (C) 2020 John W. Eaton and others. > > Octave was configured for "x86_64-w64-mingw32". > > > > >> x=linspace(-0.5,0.5,1000)' ; X=exp(1i*x.*(-10:10)) ; > > >> X(1:4,1:5) > > ans = > > > > 0.28366 - 0.95892i -0.21080 - 0.97753i -0.65364 - 0.75680i > -0.93646 - > > 0.35078i -0.98999 + 0.14112i > > 0.27405 - 0.96172i -0.21959 - 0.97559i -0.65968 - 0.75154i > -0.93889 - > > 0.34421i -0.98913 + 0.14706i > > 0.26441 - 0.96441i -0.22837 - 0.97357i -0.66568 - 0.74624i > -0.94128 - > > 0.33763i -0.98823 + 0.15300i > > 0.25474 - 0.96701i -0.23714 - 0.97148i -0.67163 - 0.74088i > -0.94362 - > > 0.33102i -0.98729 + 0.15893i > > > > >> c=X\exp(x) > > c = > > > > -0.000034734 + 0.000054543i > > 0.000612095 - 0.000907876i > > -0.005135034 + 0.007102460i > > 0.027233003 - 0.034494175i > > -0.102184601 + 0.115336800i > > 0.287266336 - 0.276216782i > > -0.622606031 + 0.467755153i > > 1.041169968 - 0.487260179i > > -1.247452051 + 0.015246341i > > 0.484987168 + 0.987560834i > > 1.412953756 - 0.175021510i > > -0.003565228 - 0.668242391i > > -0.507142432 - 0.257223380i > > 0.346131942 + 0.638760179i > > -0.148113496 - 0.545360752i > > 0.042654051 + 0.308242889i > > -0.006955253 - 0.125735396i > > -0.000093980 + 0.037054597i > > 0.000344207 - 0.007552276i > > -0.000075749 + 0.000958155i > > 0.000006065 - 0.000057233i > > > > >> norm(X*c-exp(x),"inf") > > ans = 4.6574e-15 > > >> f=@(x)exp(1i*x.*(-10:10))*c ; > > >> y=linspace(-0.5,0.5,1e6)';norm(f(y)-exp(y),"inf") > > ans = 5.2132e-15 > > >> > > > > %} > > > > ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ > > > > import numpy as np > > x=np.reshape(np.linspace(-0.5,0.5,1000),(1000,1)) > > n=np.reshape(np.arange(-10,1+10),(1,21)) > > X=np.exp((0+1j)*x.dot(n)) > > c=np.linalg.lstsq(X,np.exp(x))[0] > > np.linalg.norm(X.dot(c)-np.exp(x),ord=np.inf) > > y=np.random.rand(10**6,1)-0.5 > > np.linalg.norm(np.exp((0+1j)*y.dot(n)).dot(c)-np.exp(y),ord=np.inf) > > > > """ > > Sample run. > > Python 3.5.2 |Anaconda custom (64-bit)| (default, Jul 5 2016, 11:41:13) > > [MSC v.1900 64 bit (AMD64)] on win32 > > Type "help", "copyright", "credits" or "license" for more information. > > >>> import numpy as np > > >>> x=np.reshape(np.linspace(-0.5,0.5,1000),(1000,1)) > > >>> n=np.reshape(np.arange(-10,1+10),(1,21)) > > >>> X=np.exp((0+1j)*x.dot(n)) > > >>> c=np.linalg.lstsq(X,np.exp(x))[0] > > __main__:1: FutureWarning: `rcond` parameter will change to the default > of > > machine precision times ``max(M, N)`` where M and N are the input matrix > > dimensions. > > To use the future default and silence this warning we advise to pass > > `rcond=None`, to keep using the old, explicitly pass `rcond=-1`. > > >>> np.linalg.norm(X.dot(c)-np.exp(x),ord=np.inf) > > 3.7732211514985746e-15 > > >>> y=np.random.rand(10**6,1)-0.5 > > >>> np.linalg.norm(np.exp((0+1j)*y.dot(n)).dot(c)-np.exp(y),ord=np.inf) > > 4.815880577669688e-15 > > >>> > > """ > > > > ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ > > ---------------------------------------------------------------------- > > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
