s-watanabe-jhod / garpos-mcmc

GARPOS-MCMC is an analysis tool for GNSS-Acoustic seafloor positioning with MCMC.
GNU General Public License v3.0
6 stars 6 forks source link

Matrix not positive definite error #1

Open SquirrelKnight opened 2 years ago

SquirrelKnight commented 2 years ago

Hello, I am attempting to test this code with my data. I have quite successful results with the original GARPOS code and thought it would be worthwhile to try the MCMC edition. Here is the error I am receiving:

---------------------------------------------------------------------------
CholmodNotPositiveDefiniteError           Traceback (most recent call last)
<ipython-input-279-6c70b1f821c3> in <module>
----> 1 E_factor, E0 = E_matrix(mu_t, mu_m, ndata, mat_dt, diff_mt, same_mt, mat_tt0, fpower, rr)

/Volumes/SeaJade 2 Backup/ONC/Programs/GARPOS-MCMC/bin/garpos_mcmc_v100/eachstep.py in E_matrix(mu_t, mu_m, ndata, mat_dt, diff_mt, same_mt, mat_tt0, fpower, rr)
    100 
    101         # cholesky decomposition
--> 102         E_factor = cholesky(E0, ordering_method="natural")
    103 
    104         return E_factor, E0

sksparse/cholmod.pyx in sksparse.cholmod.cholesky()

sksparse/cholmod.pyx in sksparse.cholmod._cholesky()

sksparse/cholmod.pyx in sksparse.cholmod.Factor._cholesky_inplace()

sksparse/cholmod.pyx in sksparse.cholmod.Factor._cholesky_inplace()

sksparse/cholmod.pyx in sksparse.cholmod._error_handler()

CholmodNotPositiveDefiniteError: ../Supernodal/t_cholmod_super_numeric.c:911: matrix not positive definite (code 1)

I do not have this issue when running the first example. I have attached my SVP and observation data. Since I am not able to upload the initcfg file, I have pasted the input below.

[Obs-parameter]
Site_name = GNSSA-05
Campaign = 2021_1
Date(UTC) = 2021-08-28
Date(jday) = 2021-240
Ref.Frame = WGS84
SoundSpeed = ../../Data/GNSS-A/obs_data/GNSSA-05/GNSSA-05_2021_1-syn-svp.csv

[Data-file]
datacsv = ../../Data/GNSS-A/obs_data/GNSSA-05/GNSSA-05_2021_1-obs_sub.csv
N_shot = 67287
used_shot = 0

[Site-parameter]
Latitude0 = 48.427667
Longitude0 = -126.1744
Height0 = -21.08
Stations = GNSSA-05-01 GNSSA-05-02 GNSSA-05-03
# Array_cent :  cntpos_E    cntpos_N    cntpos_U
Center_ENU = 0.00023755002814596082 -0.02986764420476599    -403.6793396447352

[Model-parameter]
# MT_Pos :  stapos_E    stapos_N    stapos_U    sigma_E sigma_N sigma_U cov_NU  cov_UE  cov_EN
GNSSA-05-01_dPos = -155.39306204694694  370.6023875224193   -407.0126688983746  0.000   0.000   0.000   0.0 0.0 0.0
GNSSA-05-02_dPos = 399.6125768278315    -51.912264904509016 -397.01270740981886 0.000   0.000   0.000   0.0 0.0 0.0
GNSSA-05-03_dPos = -244.21880213080013  -318.7797255505246  -407.0126426260122  0.000   0.000   0.000   0.0 0.0 0.0
dCentPos = 0.0  0.0 0.0 3.0000  3.0000  3.0000  0.0 0.0 0.0
# ANT_to_TD :   forward rightward   downward    sigma_F sigma_R sigma_D cov_RD  cov_DF  cov_FR
ATDoffset = 0.0053  0.0 0.92813 0.0 0.0 0.0 0.0 0.0 0.0

Please let me know if it would be possible to help and if there is anything else that you might need from me.

GNSSA-05_2021_1-obs_sub.csv GNSSA-05_2021_1-syn-svp.csv .

SquirrelKnight commented 2 years ago

I'd like to add that I have found that by altering the rr parameter to 0.4, I am successfully able to run iterations. I am not entirely sure how this would effect the data set, however.

Additionally, it seems that because my data set is much larger, it becomes extremely slow, and unwieldy. I'm not sure if there would be a good way to address this through multiprocessing/vectorization.

s-watanabe-jhod commented 2 years ago

Thank you for using GARPOS-MCMC.

The matrix E (variance-covarinace of data error) is defined based on Eq. (14) of Watanabe et al. (2020), but multiplying a parameter rr for i /= j (for all the non-diagonal terms).

From your obs data, I guess that the "not positive definete error" comes from your GNSS-A observation system. Different from our system, your system seems to call all transponders with one signal, resulting in the same "ST" parameter for other data. It leads to the large non-diagonal term in E as defined Eq. (14) when mu_MT is close to 1. Setting smaller rr tentatively helps to suppress the non-diagonal term. That is why you succeeded in avoiding the error with small rr.

For your additional comment, I am also concerned and now developing faster scheme for the next version.

SquirrelKnight commented 2 years ago

Thank you for the response. My system does indeed send out one signal that can have multiple returns. I am using a wave glider (an autonomous GPS-A unit), which means I can survey a site for quite a long while. I am happy to hear that the data may contribute to building a faster version, and look forward to testing it.