BUG: sampler occasionally hangs

[martinjankowiak]

hello @zoj613. thanks again for the great package. i'm sometimes seeing the sampler stall/hang.

this is the smallest reproducible script i've been able to come up with. hopefully this will surface the issue on your system as well...

import numpy as np
from polyagamma import random_polyagamma

bad_state = {'bit_generator': 'PCG64', 'state': {'state': 286445528977372604878304246648472662937,
             'inc': 87136372517582989555478159403783844777}, 'has_uint32': 1, 'uinteger': 165759272}

rng = np.random.default_rng(0)
rng.__setstate__(bad_state) # if this line is commented out the sampler returns as usual

bad_z = np.array([-7.7159955279217645, -5.688383408336662, 2.4103941716465282, -4.90626055141453, -10.954245335948691, -5.773287517761037, -6.311513126139971, -7.695398984797436, -6.3832143833769655, -2.564762898556397, -2.3118010212900364, -4.718389207128771, 1.7601304947837049, -8.42732831735728, -3.862044943685949, -5.7787177292797605, -20.29932237601833, -3.6537410213424506, -8.461559085158452, -9.89944682307928, -5.696226884564253, -9.1250497748862, -7.5979075278743196, -10.757307721780466, -9.084993406805888, -7.685660133809492, -8.823505426373037, -1.165363252205105, -7.300061338287363, -6.427860810082782, -3.0221567231515736, -6.328243375567028, -1.6684474535428366, 1.2254003055223004, -6.758976991522837, -7.107719705832688, -6.981026096361312, -6.66207820518799, -6.463985934661833, -5.777396543276716, -5.18125065175635, -7.6786030941131, -7.010016457052919, -0.9000165011313692, -2.4956708005547337, -1.0192764466185116, -8.01812363871118, -3.5894788262246404, -7.522172263151698, -7.526912680584633, -2.1891046299307213, -10.424348071997604, -2.202330318247364, 2.577534253458892, -7.687467971494408, -9.688685210322241, -3.634411805182897, -1.1908965060811223, -8.582967381893358, -2.4486844411942603, -3.1510089973232405, -4.889688631223909, -8.71235438853951, 0.23921879992437667, -1.9378540362513546, -3.8802915418435457, -7.841221724612904, -0.2955492911347388, -1.3542976365653816, -6.8293614107612015, -9.210993532041117, -2.887974970403637, -7.048053453517596, -11.179085399642343, -3.201835581408629, -9.36825203899268, 1.217126729237389, -11.81015329865815, -2.8882319014328415, -3.385019075846617, -7.703097676843385, -4.220019528089502, -2.5664772197876733, -1.6384797314223714, -3.764091299623171, -3.1117359964962104, -5.210124763540389, -5.529963633124693, -12.473719197610416, -6.789538678553491, -9.939334768370019, -5.400019923735151, -11.81993965937595, -4.53756250687539, 0.5913833358729432, -0.5766533426839153, -5.843879901728336, -0.29539872879535967, -5.5296471395017255, -6.9739394740890575, -4.899325406529429, -8.190228267953618, -2.298693288379017, -12.56079570496085, -5.2361155118083476, -0.8546938483482389, -7.5999273795647815, -3.1322396929326826, -7.520184381515393, -4.382074746373018, -6.617253735425976, -13.66693607501863, -9.407616048212464, -2.4494309076866214, -2.6300391483943977, -10.863656921744937, -3.594100095197258, -0.7562807684032578, -3.0860453106408974, -8.678960150780716, -8.284999367859644, -11.270628516605488, -2.8037376000322034, -5.781402398643731, 2.697850623115234, -8.81352963529242, -0.37103120086357766, -5.722245592378262, -4.406450954903647, -9.880914623244749, -2.0592691145346746, -5.058151894855059, 0.9438105597141462, -2.013557031693369, -2.2185189585681138, -7.9966281089095, 3.489025571184194, -3.279924879114724, -1.186609088971064, -9.970855880293058, -6.667642513611053, -8.84694390533457, -3.978299814282288, -6.06052387224828, -3.0605332422625096, -5.049904760405894, -7.05815944833525, -2.2529170769032745, -10.89285943942715, -11.068038514304959, -2.3860733591154197, -0.08798811481825997, -9.316088255444871, -6.4668465057920965, -6.811544859596736, -4.05258596226946, -5.280669655258182, -7.365182715882279, -4.543439696164741, -5.148616554927539, -3.750774897823457, -12.193567338305948, -0.6957150084697634, 2.808564255056562, -0.6876715432885856, -8.02318701562676, -6.977896662944611, -9.524408581426663, -9.781958289842166, -4.271102793111396, -1.0411750040334482, -8.714375164603663, -5.21215745145907, -3.963140641101465, -2.093652102829111, -8.136657864869417, -3.7372133579288045, -4.544260411064354, -2.315289889567245, -10.575288150277192, -10.836499727310812, -2.2017240079309413, -4.319009311151649, -3.393951127373044, 2.3319051665861164, -6.415511447193612, -12.000559977123228, -3.932846640222944, -4.2839988799104525, 0.3748367571387474, -4.38119258087972, -4.9975878200024155, -1.215029777311032, -4.432040728644365, -4.055317480451655, -8.015277250238704, -7.791193041364883, -1.208604812146648, -7.464106652678793, -3.3089459574695512, 0.565354111599949, 4.006234308066209, -5.411213095562011, -5.333275728492514, -7.168529856598479, -4.888233382212302, -11.264199152466603, -6.165752583907858, -7.63420257164309, -3.0591368035023128, -4.1187622907820485, -15.213779959499778, -6.1490315083215155, -5.3856648070059325, -2.908478940128666, -11.145910344518146, -6.997522954584618, -7.299823673363338, -8.129460491657383, -0.9268450449842698, -4.477811838108732, -7.464161846723406, -6.660454515563282, 0.08286648798164276, -3.0350911859177643, -5.231581623118869, -4.0344359736463336, -3.007616728277827, -6.812590583179762, -4.751400902831156, -6.2289966072471685, -6.698795132832843, -6.046464323967813, -9.166306216389106, -6.749800617887599, 1.0684002553397738, 3.035042848011444, -8.956817926179387, -5.432044283421149, -2.828781674885361, -3.0467644287247553, 1.1394259540556186, 1.7625068559757242, -10.022373510387043, -6.256482949081717, -3.641392297856366, -6.971862587394393, -3.5125729470562996, -9.911990130316614, -0.37314856765888305, -6.679089485721889, -3.562972669877401, -1.6164110624843309, -6.353916245342262, -3.572645464566901, -8.919922783159063, -6.1019248023136585, -0.1555105784862203, -1.828608864388647, -4.177580353542247, -5.612723250588665, 1.2993931058958985, -9.570976847117686, -3.8136562734762087, -2.8270100999665284, -5.742013130680353, -9.883000701084573, -7.32928907338796, -2.944968201160664, -5.929924037218706, -3.7345046817224348, -4.22621858869544, -8.378613350426871, 0.7021251460062743, -4.100515471133535, -11.943186161950807, -8.339145960195191, 1.6986542094024006, -7.316982269280614, -7.4089735116666064, -3.0046255099671892, -2.5873212081597377, -6.827621059495563, -1.438531959738429, 1.8796104554474073, 1.5050618611081399, -9.028777963657493, -0.6891704018741098, -6.837056055237005, -5.583166831901322, -0.08426140515610658, -2.8892378920633335, -2.3993629823340865, -6.763798068451533, 3.288598191811139, -8.206104540496137, -4.145173066015854, -4.227512041516999, -10.535386674255394, -8.16948303108025, -5.372696166836505, 0.6429787585704672, -8.122691482350891, -0.30523728485417667, -7.237694943445989, -0.13473236220564644, -5.437492787011896, -2.564590189912019, -5.396626372177929, -7.1683511306506436, 1.1849062433080482, -7.92707764316654, -11.0703241395689, -6.555616270263718, -9.497107742926442, -6.2114905714677455, -7.28886644713213, -6.914132906415123, -6.723028005985412, -6.474677110329841, -2.533929317578626, -7.525411506029455, -1.3099847524292487, -7.554562569783455, -0.9813638358166576, -7.751284491848235, -8.817102676207508, -9.370984076524948, -5.6011046151125345, -2.5060527557251966, -2.3914963973119985, -12.646179375261301, -4.919562103924632, -1.2997636152623437, -5.308285827321867, -2.359816554352985, -2.8208051952601227, -3.12147881303965, -2.5648873878065426, -4.747637220537003, -2.624172953883202, -9.950726813032155, -7.141427173562353, -6.719985768830318, -9.472346941311162, -11.2633340949679, -11.262049026248997, -1.0490965770517695])

random_polyagamma(np.ones(bad_z.shape), bad_z, random_state=rng)

[zoj613]

Hi @martinjankowiak thank you once again for reporting any bugs you come across. This is a great opportunity for me to fix any weaknesses in the implementation that I was not aware of. Since I have no idea how many people are using this package I only have had my own use cases as a way of "testing" the library's quality.

I can definitely reproduce this issue locally, and it sure does go away when I remove the __setstate__ line. Based on the bad state used, the only method affected by it is the devroye algorithm. The other 2 return just fine without hanging. So I did a read-through of the module to see what could cause it to "hang". Then it finally hit me: I had removed a crucial "is-positive" check for generated proposals in a previous commit (d014bae#diff-3ce312a6654d7849f3c0ab397696d5f716cce2505e05ba24833d22ce08e2eb10L31-L37) because it seemed useless at face-value and that I was being too careful. See:

polyagamma/src/pgm_devroye.c

Lines 27 to 46 in da31ac8

	/*
	* Compute a_n(x|t), the nth term of the alternating sum S_n(x|t)
	*
	* NOTE
	* ----
	* pr->x is guaranteed to be always positive due to proposal support so no
	* need for extra checks other than to test if its greater than the truncation
	* point.
	*/
	static PGM_INLINE float
	piecewise_coef(int n, parameter_t const* pr)
	{
	if (pr->x > T) {
	double b = PGM_PI * (n + 0.5);
	return (float)b * expf(-0.5 * pr->x * b * b);
	}
	double a = n + 0.5;
	return expf(-1.5 * (PGM_LOGPI_2 + pr->logx) - 2. * a * a / pr->x) *
	(float)(PGM_PI * a);
	}

I see that a very small -x or x=0 can be generated as a proposal sample during the accept-rejection sampling which can throw off the sampler and cause it to hang. So I re-introduced the check for positive x in the piecewise_coef function and then voila! The hang disappeared.

Because re-introduction of that line should not affect the sequence of the states, then I think we can be confident that the cause of the hang was likely that a value of x=0 was generated and not handled accordingly during the sampling.

[zoj613]

A little extra digging shows the source of this issue with the proposal sample x being outside the positive range of numbers. The culprit is this line:

polyagamma/src/pgm_devroye.c

Line 133 in da31ac8

x = w - sqrt(w * w - 1 / pr->z2);

When w * w and 1 / pr->z2 are very close, the difference can have a negative sign, which results in x being -nan because the sqrt of a negative number does not exist.
e.g:

w*w=0.0671855778511380, 1/z2=0.0671855778511380, diff=-0.0000000000000000

The above is the standard output of the values in the line where it happens. According to my findings, it happens for the element z=-6.719985768830318 in your input for the given initial state, which is the 5th element from the end in the array bad_z (but this info isn't relevant since the input itself is not the issue but the way implementation does not account for the rare case of having these two terms w*w and 1/z2 being equal).

So it seems like real solution is to ensure that the sign of the expression inside the sqrt is not negative in the case that it evaluates to zero.

A bit of algebra shows that the expression inside the sqrt can be re-written as:
y*y * (z + 0.25 * y*y) / (z2 * z2), which wont suffer from the issues above due to taking a difference of equal terms. The expression simplifies to: (y/z)^2 * (1/z + 0.25 * (y/z)^2)

[zoj613]

@martinjankowiak I've put out an alpha release for your immediate use at https://pypi.org/project/polyagamma/1.3.2a0/. Let me know if you run into any issues. Hopefully #80 can be solved soon so I can put out an official release.