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
