%. is for _linear_ least squares, so you have to convert a non-linear model
to a linear one, if you can. (If you can't, %. won't help you.) Your
model is:
y = a * ^ b*x
Therefore:
(^.y) = (^.a) + (^. ^ b*x)
(^.y) = (^.a) + (b*x)
So for the data to be fit, compute c=: (^.y) %. 1,.x , whence a=^{.c and
b={:c.
On Wed, Apr 28, 2021 at 9:17 PM 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
