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
