Least square fitting by SVD

svd fitting

cc by Shigeto R. Nishitani, 2023

簡単な最小二乗法によるfittingでも計算誤差が出てしまうという例です．0からだいぶ離れたところで，微妙なエネルギー差をfittingしようとして出てきました．

結論はSVD(single value decomposition, 特異値分解)を使えというNumRecipeの教えを再確認．．．やれやれ．


# 最初のうまくいってない結果

-   include NumericCalcs
-   fit
-   plot


## Ruby module for Numerical calculations

```ruby
  require "matrix"
  require "rbplotly"

  module NumericCalcs
    include Math

    def fitting(x, d_y, dim = 3)
      y = d_y

      # make design matrix
      n, m = x.size, dim

      av = Matrix.zero(n, m)
      yy = Vector.zero(n)
      n.times do |i|
        m.times do |j|
          av[i,j] = x[i] ** j
        end
        yy[i] = y[i]
      end

      # calc inverse for non-square matrix
      ai = (av.transpose * av).inv
      b = av.transpose * yy
      return ai * b
    end

    def polynomials(coeffs, x)
      sum = 0.0
      coeffs.each_with_index do |c, i|
        sum += c * x ** i
      end
      sum
    end

    # x_ticks
    def mk_x(init, fin, dev) 
      x = []
      (init..fin).step(dev) { |i| x << i }
      x
    end

    def plot_data_and_func(xx, y_vals, func)
      dev = (xx.max - xx.min) / 100.0
      xx2 = mk_x(xx.min, xx.max, dev)

      traces = [{ x: [], y: [], name: "data" },
                { x: [], y: [], name: "fitted func" }]

      xx.size.times do |i| # rev for Numo::DFloat
        traces[0][:x] << xx[i]
        traces[0][:y] << y_vals[i]
      end

      xx2.each do |x|
        traces[1][:x] << x
        traces[1][:y] << func.call(x)
      end
      pl = Plotly::Plot.new(data: traces)
      Plotly.default_auth()
      return pl
    end
  end

  include NumericCalcs
```

    #<Module:0x000000013e3cd518>

```ruby
  p x_vals = [950, 960, 970, 980, 990]
  p y_vals = [-4.177412829697432, -4.1848393364765, -4.192296922237924, 
      -4.199785335389045, -4.207304330283102]
```

    [950, 960, 970, 980, 990]
    [-4.177412829697432, -4.1848393364765, -4.192296922237924, -4.199785335389045, -4.207304330283102]
    
    
    
    
    
    [-4.177412829697432, -4.1848393364765, -4.192296922237924, -4.199785335389045, -4.207304330283102]


## main results

```ruby
  xx=x_vals
  yy = y_vals
  p vals=fitting(xx, yy, 4)

  f_plot = lambda { |x| polynomials(vals, x) }
  pl = plot_data_and_func(xx, yy, f_plot)
  pl.show
  #pl.savefig('./figs/fit_plot0.png')
```

    Vector[-3.650146484375, -0.0003299713134765625, -2.7194619178771973e-07, 4.1382008930668235e-11]

<a id="orgb39bd82"></a>

![./figs/fit\_plot0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151211/05b268ec-f35d-debb-bbac-9b33d3ab6f35.png)

というのが default 計算の結果． これを改善できないかという話．


# 少し試行錯誤

明らかにずれていて，Maple の fitting を使うとばっちり．NumericCalcs::fitting が甘いのはあきらかだが，どうすれば．．．

```ruby
  av = Matrix[[1, 950, 902500, 857375000], [1, 960, 921600, 884736000], [1, 970, 940900, 912673000], [1, 980, 960400, 941192000], [1, 990, 980100, 970299000]]
  yy = Vector[-4.177412829697432, -4.1848393364765, -4.192296922237924, -4.199785335389045, -4.207304330283102]

  ai = (av.transpose * av).inv
  b = av.transpose * yy
  ai * b
```

    Vector[-3.650146484375, -0.0003299713134765625, -2.7194619178771973e-07, 4.1382008930668235e-11]

すると確かに．おなじ．

でも，これを Maple でやるとまた違った答え

    [-3.650634765625, -0.0003314018249511719, -0.2756714821*10^(-6), 0.4160938261*10^(-10)])

になる．これは，数値誤差(Ruby の inv かな)のようなんで，放置．


## 1000->1 に(効果なし)

```ruby
  x_vals2 = [950, 960, 970, 980, 990]
  x_vals2.map!{|val| val/1000.0}
  y_vals = [-4.177412829697432, -4.1848393364765, -4.192296922237924, 
      -4.199785335389045, -4.207304330283102]
```

    [-4.177412829697432, -4.1848393364765, -4.192296922237924, -4.199785335389045, -4.207304330283102]

```ruby
  xx=x_vals2
  yy = y_vals
  vals=fitting(xx, yy, 4)

  f_plot = lambda { |x| polynomials(vals, x) }
  pl = plot_data_and_func(xx, yy, f_plot)
  pl.show
  #pl.savefig('./figs/fit_fail1.png')
```

<a id="orge5e9f59"></a>

![./figs/fit\_plot1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151211/2d824ffa-6533-eae8-36ae-74a026a47dd8.png)

ちょっと変わりますが，差は大きいです．


# Ruby の svd


## Numo の導入

steep なんやろな．なら，SVD だが．Ruby にあったっけ？

```sh
  > gem install ruby-svd
  > emacs /Users/bob/.rbenv/versions/3.0.2/lib/ruby/gems/3.0.0/extensions/arm64-darwin-20/3.0.0/ruby-svd-0.5.1/gem_make.out
```

うーーん．gccのcompile エラー．

-   [numo-linalgを mac で動かしてみた](https://qiita.com/daddygongon/items/73a83565b7969ffc5d1a)

numo がそのまま動いた．へぇ，やるね．2年ほど前に動かなんでNumoを使ったコードをMatrixに書き換えたんですが．．．

どこが対応したんやろ．．． -[Adddynamic\_lookup linker option for macOS and Ruby3.1](https://github.com/ruby-numo/numo-linalg/blob/master/ext/numo/linalg/lapack/extconf.rb)と４ヶ月ほど前に修正したはります

```ruby
  require 'numo/gnuplot'
  require 'numo/narray'
  require 'numo/linalg'

  x, y=x_vals, y_vals

  # make design matrix
  n = 5 # data
  m = 4 # dimension
  av = Numo::DFloat.zeros(n, m)
  yy = Numo::DFloat.zeros(n)
  n.times do |i|
    m.times do |j|
      av[i,j]=x[i]**j
    end
    yy[i] = y[i]
  end

  # fit by using matrix inv
  p ai = Numo::Linalg.inv(Numo::Linalg.dot(av.transpose,av))
  p b = Numo::Linalg.dot(av.transpose,yy)
  $vars = Numo::Linalg.dot(ai,b)
  p ['$vars', $vars]

  # block for lazy evaluation in plot_data_and_func
  free_e_plot = lambda { |x| polynomials($vars, x) }

  pl = plot_data_and_func(x, y, free_e_plot)
  pl.show
```

    Numo::DFloat#shape=[4,4]
    [[5.78259e+10, -1.7888e+08, 184428, -63.3749], 
     [-1.7888e+08, 553356, -570.525, 0.196052], 
     [184428, -570.525, 0.588234, -0.00020214], 
     [-63.3749, 0.196052, -0.00020214, 6.94638e-08]]
    Numo::DFloat#shape=[4]
    [-20.9616, -20333.5, -1.97284e+07, -1.91454e+10]
    ["$vars", Numo::DFloat#shape=[4]
    [-3.64941, -0.000331879, -2.73809e-07, 4.18368e-11]]

<a id="orgc66f34b"></a>

![./figs/fit\_plot2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151211/d95b0c33-c8e3-2252-a2cb-5e4c8465c875.png)

少し変わってるのも(Matrix計算精度が)気になるけど，まあ，変わらんはね．


## Numoの動作確認

SVDでの実装ですが，先ずは，[講義のサイト](https://nbviewer.org/github/daddygongon/jupyter_num_calc/blob/master/numerical_calc/least_square_fit.ipynb#特異値分解(Singular-Value-Decomposition)による解)から確実に動くのを取り出して，確認．

```ruby
  x = Numo::NArray[[1,2,3,4]]
  y = Numo::NArray[[0,5,15,24]]

  fitting(x,y)
```

    Vector[(-9/2), (16/5), (1/1)]

```ruby

  x = Numo::NArray[[1,2,3,4]]
  y = Numo::NArray[[0,5,15,24]]

  # make design matrix
  n, m = 4, 3 # data, dimension
  av = Numo::DFloat.zeros(n, m)
  yy = Numo::DFloat.zeros(n)
  n.times do |i|
    m.times do |j|
      av[i,j]=x[i]**j
    end
    yy[i] = y[i]
  end

  # fit by using matrix inv
  require 'numo/linalg'
  ai = Numo::Linalg.inv(Numo::Linalg.dot(av.transpose,av))
  b = Numo::Linalg.dot(av.transpose,yy)
  vars = Numo::Linalg.dot(ai,b)
  p vars

  # block for lazy evaluation in plot_data_and_func
  y_fit = lambda { |x| polynomials(vars, x) }
  pl = plot_data_and_func(x, yy, y_fit)
  pl.show
```

    Numo::DFloat#shape=[3]
    [-4.5, 3.2, 1]

<a id="org21dc92d"></a>

![./figs/fit\_plot3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151211/6917049f-c7a0-a654-73a5-52ac63f6e8af.png)


## Numo の svd の出力順は python と違う

[numpy.linalg.svd](https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html)にある通り，

    a : (m,n array)
    u,s,vh = linalg.svd(a)で
    u : (m , m)
    s : (k) よくわかってないよね．．．
    vh : (n, n)

ところが Numo::Linalg.svd では，

    s,u,vh = Numo::LInalg.svd(a)
    s : k
    u : m x m
    vh : n x n

という感じ．

これに関するマニュアルはなさそう．[ソース](https://github.com/ruby-numo/numo-linalg/blob/master/lib/numo/linalg/function.rb)が

```ruby
  l.323 Lapack.call(:gesdd, a, jobz:job)[0..2]
```

となってなるんで，s,u,vの順で正しそう．vは$v^T$とか書いてあったような．．． 確認して使わないと．

```ruby
  ss, uu, vs =  Numo::Linalg.svd(av)

  iS = Numo::DFloat.zeros(m, n)
  m.times do |j|
      iS[j,j]=1.0/ss[j]
  end
  p iS

  p left = Numo::Linalg.dot(vs.transpose,iS)
  p right = Numo::Linalg.dot(uu.transpose,yy)

  Numo::Linalg.dot(left, right)
```

    Numo::DFloat#shape=[3,4]
    [[0.0509649, 0, 0, 0], 
     [0, 0.584088, 0, 0], 
     [0, 0, 3.75583, 0]]
    Numo::DFloat#shape=[3,4]
    [[-0.00421136, -0.393889, 2.75587, 0], 
     [-0.0138812, -0.404752, -2.50719, 0], 
     [-0.0488569, 0.14895, 0.474791, 0]]
    Numo::DFloat#shape=[4]
    [-28.6152, 1.83564, -1.41424, 1.34164]
    
    
    
    
    
    Numo::DFloat#shape=[3]
    [-4.5, 3.2, 1]

係数は，invを使ったのからsvdを使ったのに変えても一致．


## 実働の確認

元のデータでinvを使ったのをまずは確認．

```ruby
  x, y=x_vals, y_vals

  # make design matrix
  n = 5 # data
  m = 4 # dimension
  av = Numo::DFloat.zeros(n, m)
  yy = Numo::DFloat.zeros(n)
  n.times do |i|
    m.times do |j|
      av[i,j]=x[i]**j
    end
    yy[i] = y[i]
  end

  # fit by using matrix inv
  require 'numo/linalg'
  ai = Numo::Linalg.inv(Numo::Linalg.dot(av.transpose,av))
  b = Numo::Linalg.dot(av.transpose,yy)
  $vars = Numo::Linalg.dot(ai,b)
  p ['$vars', $vars]

  # block for lazy evaluation in plot_data_and_func
  free_e_plot = lambda { |x| polynomials($vars, x) }

  pl = plot_data_and_func(x, y, free_e_plot)
  pl.show
```

    ["$vars", Numo::DFloat#shape=[4]
    [-3.64941, -0.000331879, -2.73809e-07, 4.18368e-11]]

<a id="orgd0f0f59"></a>

![./figs/fit\_plot2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151211/d95b0c33-c8e3-2252-a2cb-5e4c8465c875.png)

```ruby
  ss, uu, vs =  Numo::Linalg.svd(av)

  iS = Numo::DFloat.zeros(m, n)
  m.times do |j|
      iS[j,j]=1.0/ss[j]
  end
  #p iS
  left = Numo::Linalg.dot(vs.transpose,iS)
  right = Numo::Linalg.dot(uu.transpose,yy)

  p vars = Numo::Linalg.dot(left, right)

  # block for lazy evaluation in plot_data_and_func
  free_e_plot = lambda { |x| polynomials(vars, x) }

  pl = plot_data_and_func(x, y, free_e_plot)
  pl.show
```

    Numo::DFloat#shape=[4]
    [-3.65027, -0.000331299, -2.74726e-07, 4.14365e-11]

<a id="orgfaae2ae"></a>

![./figs/fit\_plot4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151211/0ec732f1-aafc-3537-3c4f-e1d0f6320bab.png)

さすが svd.


# 最終形（完成版fitting\_svd)

```ruby
  require 'numo/gnuplot'
  require 'numo/narray'
  require 'numo/linalg'

  def fitting_svd(x, d_y, dim = 3)
    y = d_y

    # make design matrix
    n, m = x.size, dim
    av = Numo::DFloat.zeros(n, m)
    yy = Numo::DFloat.zeros(n)
    n.times do |i|
      m.times do |j|
        av[i,j]=x[i]**j
      end
      yy[i] = y[i]
    end
    ss, uu, vs =  Numo::Linalg.svd(av)

    iS = Numo::DFloat.zeros(m, n)
    m.times do |j|
      iS[j,j]=1.0/ss[j]
    end
    #p iS
    left = Numo::Linalg.dot(vs.transpose,iS)
    right = Numo::Linalg.dot(uu.transpose,yy)

    return Numo::Linalg.dot(left, right)
  end
```

    :fitting_svd

```ruby
  p x_vals = [950, 960, 970, 980, 990]
  p y_vals = [-4.177412829697432, -4.1848393364765, -4.192296922237924, 
      -4.199785335389045, -4.207304330283102]
      
  fitting_svd(x_vals, y_vals, 4)
```

    [950, 960, 970, 980, 990]
    [-4.177412829697432, -4.1848393364765, -4.192296922237924, -4.199785335389045, -4.207304330283102]
    
    
    
    
    
    Numo::DFloat#shape=[4]
    [-3.65027, -0.000331299, -2.74726e-07, 4.14365e-11]

```ruby
  vars = fitting_svd(x_vals, y_vals, 4)
  free_e_plot = lambda { |x| polynomials(vars, x) }

  pl = plot_data_and_func(x_vals, y_vals, free_e_plot)
  pl.show
```

<a id="org6f5496e"></a>

![./figs/fit\_plot4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151211/0ec732f1-aafc-3537-3c4f-e1d0f6320bab.png)


# Mapleとの比較(overtaken?!)

得られた結果をMapleの結果と比較すると，

```octave
  f1 := x -> 7.526568118*x^3/10^9 + 
  (-1)*0.00002204463465*x^2 + 
  0.02077149587*x - 
  10.46814903

  f2:=unapply(-3.65027 
  -0.000331299*x 
  -2.74726e-07*x^2
  + 4.14365e-11*x^3, x);

  eV2J:=(1.60218*10**(-19))*(6.0221367*10**(23));
  subs(x=970,-x*diff(f1(x)*eV2J, x,x)); # by Maple fitting
  subs(x=970,-x*diff(f2(x)*eV2J, x,x)); # by Numo::svd
                        eV2J := 96485.46978
                            26.63998203
                            28.85328724
```

とズレがある．これをより広域のplotで比べると

![./maple\_figs/free\_energy\_fit\_compare.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/151211/7b1a0ade-49bc-678e-3c0c-910773933a19.png)

Fitted free energy by SVD(blue) vs Maple(red).

うーーーん．見た目は Numo::Linalg.svdの勝ちだね．

実際，少しずれた温度(i.e. $T$=980) で$F(T)$を比べるとMapleは全然違う値．

エントロピー($S$)は $$Z = \left(\frac{\exp(-\Theta_E/2T)}{1-\exp(-\Theta_E/T)} \right)^{3N} \\S = k_B \ln Z$$ だな．


# Qiita記事作成のメモ

fit\_plot.ipynbより途中で，[qiita\_org](https://qiita.com/yamatoken/items/9314161758e778fd82c0)

```sh
  > qiita upload fit_plot.org
```

でfigsを登録したのちに，以下でQiita記事を作成．

```sh
  > jupyter nbconvert --to markdown fit_plot.ipynb
  > pandoc -f markdown -t org -o fit_plot.org fit_plot.md
  > rake fit_plot > tmp.org
  > qiita post tmp.org
```



------
 - **source** ~/current_research/higher_moment/moment_al/einstein/harmonic_free/tmp.org