Python numpy : linalg.pinv() 太不精确
Python numpy : linalg.pinv() too imprecise
我最近一直在算法中使用 numpy 矩阵，但遇到了一个问题：
我一共用了3个矩阵
m1 = [[ 3 2 2 ... 3 2 3]
[ 3 3 3 ... 2 2 2]
[500 501 502 ... 625 626 627]
...
[623 624 625 ... 748 749 750]
[624 625 626 ... 749 750 751]
[625 626 627 ... 750 751 752]]
m1 是 (128,128) 奇异 方阵。前两行看似随机的 2 和 3 序列。接下来的行按照算法从 500 开始填充，从第三行第一列开始为每一行和每一列加一个。
m2 = [[ 2 3 500 ... 623 624 625]
[ 2 2 500 ... 623 624 625]
[ 3 2 500 ... 623 624 625]
...
[ 2 3 500 ... 623 624 625]
[ 2 2 500 ... 623 624 625]
[ 3 2 500 ... 623 624 625]]
m2也是一个(128,128)奇异方阵。这次，随机序列归因于前两列。每行的其余部分填充500、501、502、503等。
m3 = [[ 790 784 157500 ... 196245 196560 196875]
[ 804 811 161000 ... 200606 200928 201250]
[ 180501 180411 36064000 ... 44935744 45007872 45080000]
...
[ 219861 219771 43936000 ... 54744256 54832128 54920000]
[ 220181 220091 44000000 ... 54824000 54912000 55000000]
[ 220501 220411 44064000 ... 54903744 54991872 55080000]]
m3 = m1*m2
所以我想做的是使用 m1 和 m3 恢复 m2。理论上，我所要做的就是执行下面的代码m2 = (m1**-1)*m3。不幸的是，由于 m1 是奇异矩阵，无法计算其逆矩阵，即使可以，矩阵也太大，导致许多数值不精确。
相反，我决定使用 m1 的 Moore-Penrose Inverse，它不需要矩阵是非奇异的，并且与逆矩阵类似，理论上可以恢复 m2, 使用 np.linalg.pinv(m1) * m3.
我再次使用“理论上”一词，因为事实证明，numpy 在涉及大矩阵的此类计算时过于不精确，这是我获得的结果 m2:
[[ 2.46207616 2.48959603 500. ... 623. 624.
625. ]
[ 2.38612549 2.61197086 500. ... 623. 624.
625. ]
[ 2.38711085 2.6125801 500. ... 623. 624.
625. ]
...
[ 2.61998539 2.37184747 500. ... 623. 624.
625. ]
[ 2.54403472 2.4942223 500. ... 623. 624.
625. ]
[ 2.62195611 2.37306595 500. ... 623. 624.
625. ]]
如您所见，m1 的整个 'filler' 部分计算正确，没有问题。但是，前两列似乎有问题，将数字四舍五入为 2s 和 3s 给我一个不正确的 m2.
我正在寻找一种方法，使 np.linalg.pinv() 方法的浮点计算更加精确，以便它可以获得序列的正确值，因为这些值非常重要。
通过一些研究，我了解到 np.linalg.pinv() 有一个名为 rcond 的参数，描述如下：
rcond : (…) array_like of float
Cutoff for small singular values. Singular values smaller (in modulus) than rcond * largest_singular_value (again, in modulus) are set to zero. Broadcasts against the stack of matrices
rcond 默认设置为 1e-15。我认为进一步减少这个数字可能有助于不精确。 1e-16 还不够，从 1e-17 开始，我得到了非常奇怪的值，例如：
[[ 3.000e+00 3.000e+00 5.000e+02 ... 6.230e+02 6.240e+02 6.250e+02]
[ 1.100e+01 4.000e+00 1.840e+02 ... 1.722e+03 2.032e+03 1.831e+03]
[-3.000e+00 -5.000e+00 -4.030e+02 ... -1.232e+03 -7.400e+02 -1.272e+03]
...
[ 1.100e+01 1.200e+01 2.164e+03 ... 4.030e+03 4.872e+03 1.873e+03]
[-1.200e+01 -9.000e+00 -1.618e+03 ... -3.240e+03 -2.519e+03 -4.167e+03]
[ 2.000e+01 2.600e+01 4.535e+03 ... 5.165e+03 5.881e+03 5.189e+03]]
所以，基本上，我被困住了，我不知道如何提高精度。
最糟糕的是，我有一个可以显着提高浮点精度的模块，它称为 mpmath 并且还有一些矩阵似乎更适合我的算法，如 numpy 矩阵。但是 mpmath 没有计算伪逆的方法，并且 numpy 没有将自己的浮点精度调整为 mpmath.
设置的值
你有什么建议我可以尝试使用伪逆方法得到正确的 m2 吗？
你的烦恼与pinv准确与否无关
正如您自己注意到的那样，您的矩阵严重缺乏秩，m1 的秩为 4 或更低，m2 的秩为 3 或更低。因此，您的系统 m1@x = m3 极度欠定，无法恢复 m2.
即使我们把我们知道的关于 m2 结构的所有信息都扔进去，即前两列 3 和 2，其余 500 向上计数，仍然有大量的组合解。
如果有足够的时间，下面的脚本可以找到它们。在实践中，我没有看到超过 32x32 矩阵，在下面显示的 运行 中产生
15093381006 个不同的有效重建 m2' 满足 m1@m2' = m3 和我刚才提到的结构约束。
import numpy as np
import itertools
def make(n, offset=500):
offset -= 2
res1 = np.add.outer(np.arange(n), np.arange(offset, offset+n))
res1[:2] = np.random.randint(2, 4, (2, n))
res2 = np.add.outer(np.zeros(n, int), np.arange(offset, offset+n))
res2[:, :2] = np.random.randint(2, 4, (n, 2))
return res1, res2
def subsets(L, n, mn, mx, prepend=[]):
if n == 0:
if mx >= mn:
yield prepend
elif n == 1:
for l in L[L.searchsorted(mn):L.searchsorted(mx, 'right')]:
yield prepend + [l]
else:
ps = L.cumsum()
ps[n:] -= ps[:-n]
ps = ps[n-1:]
for i in range(L.searchsorted(mn - np.sum(L[len(L)-n+1:])),
ps.searchsorted(mx, 'right')):
yield from subsets(L[i+1:], n-1, mn - L[i], mx - L[i],
prepend = prepend + [L[i]])
def solve(m1, m3, ci=0, offset=500):
n, n = m1.shape
col = m3.T[ci]
n3s = col[3] - col[2] - 2 * n
six = col[2] - offset * (col[3] - col[2]) - n * (n-1)
idx = np.lexsort(m1[:2])
m1s = m1[:2, idx]
sm = m1s[1].searchsorted(2.5)
sl = m1s[0, :sm].searchsorted(2.5)
sr = sm + m1s[0, sm:].searchsorted(2.5)
n30 = n - sl - sr + sm
n31 = n - sm
n330 = col[0] - 4*n - 2*n3s - 2*n30
n331 = col[1] - 4*n - 2*n3s - 2*n31
for n333 in range(max(0, n330 - sm + sl, n331 - sr + sm),
min(n - sr, n330, n331) + 1):
n332 = n330 - n333
n323 = n331 - n333
n322 = n3s - n332 - n323 - n333
mx333 = six - idx[sl:sl+n332].sum() - idx[sm:sm+n323].sum() \
- idx[:n322].sum()
mn333 = six - idx[sm-n332:sm].sum() - idx[sr-n323:sr].sum() \
- idx[sl-n322:sl].sum()
for L333 in subsets(idx[sr:], n333, mn333, mx333):
mx332 = six - np.sum(L333) - idx[sm:sm+n323].sum() \
- idx[:n322].sum()
mn332 = six - np.sum(L333) - idx[sr-n323:sr].sum() \
- idx[sl-n322:sl].sum()
for L332 in subsets(idx[sl:sm], n332, mn332, mx332,
prepend=L333):
mx323 = six - np.sum(L332) - idx[:n322].sum()
mn323 = six - np.sum(L332) - idx[sl-n322:sl].sum()
for L323 in subsets(idx[sm:sr], n323, mn323, mx323,
prepend=L332):
ex322 = six - np.sum(L323)
yield from subsets(idx[:sl], n322, ex322, ex322,
prepend=L323)
def recon(m1, m3, ci=0, offset=500):
n, n = m1.shape
nsol = nfp = 0
REC = []
for i3s in solve(m1, m3, ci, offset):
rec = np.full(n, 2)
rec[i3s] = 3
if not np.all(m3.T[ci] == m1@rec):
print('!', rec, m3.T[ci], m1@rec)
nfp += 1
else:
nsol += 1
REC.append(rec)
print('col', ci, ':', nsol, 'solutions,', nfp, 'false positives')
return np.array(REC)
def full_recon(m1, m3, offset=500, subsample=10):
n, n = m1.shape
col0, col1 = (recon(m1, m3, i, offset) for i in (0, 1))
yield col0.shape[0], col1.shape[0]
if not subsample is None:
col0, col1 = (col[np.random.choice(col.shape[0], subsample)]
if col.shape[0] > subsample else col
for col in (col0, col1))
print('col 0', col0)
print('col 1', col1)
for c0, c1 in itertools.product(col0, col1):
out = np.add.outer(np.zeros(n, int), np.arange(offset-2, offset+n-2))
out[:, :2] = np.c_[c0, c1]
yield out
def check(m1, m3, offset=500, subsample=10):
for cnt, m2recon in enumerate(full_recon(m1, m3, offset, subsample)):
if cnt == 0:
tot0, tot1 = m2recon
continue
assert np.all(m3 == m1@m2recon)
print(cnt, 'solutions verified out of', tot0, 'x', tot1, '=', tot0 * tot1)
样本运行:
>>> m1, m2 = make(32)
>>> check(m1, m1@m2)
col 0 : 133342 solutions, 0 false positives
col 1 : 113193 solutions, 0 false positives
col 0 [[2 2 3 2 2 2 3 2 3 3 3 3 3 2 3 3 2 2 2 2 2 3 3 3 2 2 2 2 2 2 2 2]
[2 2 3 2 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 3 3 3 2 2 2 3 2 3 2 2 2 2]
[2 3 3 2 3 3 2 2 2 3 2 3 3 2 2 2 2 3 3 3 2 2 2 2 3 2 2 2 2 2 2 3]
[2 2 3 3 2 2 3 3 3 2 2 3 2 3 3 3 2 2 2 2 2 2 2 3 2 3 3 2 2 2 2 2]
[2 3 3 3 2 3 3 2 2 2 2 3 2 2 3 3 2 2 3 2 2 3 2 2 2 2 2 2 3 3 2 2]
[3 2 2 2 3 2 3 3 3 2 3 2 2 2 2 3 3 3 2 2 2 2 3 3 2 3 2 2 2 2 2 2]
[2 2 2 3 3 3 2 2 3 3 3 3 3 2 2 2 2 3 2 2 2 3 2 2 2 3 2 2 3 2 2 2]
[2 2 2 3 3 2 3 3 3 2 3 2 2 3 3 3 2 3 2 2 2 2 2 2 2 2 2 3 2 3 2 2]
[3 2 3 3 2 2 3 2 2 3 2 3 3 2 2 2 3 2 3 2 2 2 2 3 2 3 2 2 3 2 2 2]
[3 2 2 3 3 2 3 3 2 2 2 3 2 3 2 2 3 2 3 3 2 2 2 2 2 3 2 2 2 2 2 3]]
col 1 [[2 2 2 2 3 3 3 3 3 3 3 3 2 2 3 2 2 2 3 3 2 3 2 2 2 3 2 3 2 2 2 3]
[2 3 3 2 3 3 2 2 2 3 3 3 2 2 3 3 2 3 2 3 2 2 3 2 2 2 2 3 3 2 2 3]
[3 2 3 2 2 3 3 2 2 3 3 3 2 2 3 2 3 2 3 2 3 2 3 3 2 2 2 2 2 3 3 2]
[2 2 3 2 3 3 3 2 3 3 2 3 2 3 2 3 2 3 2 2 3 2 3 2 2 3 2 3 2 2 2 3]
[3 3 2 2 2 2 3 2 3 3 3 3 2 3 3 2 2 3 3 2 2 2 2 2 3 2 2 3 3 3 2 2]
[3 3 2 3 2 2 2 3 3 3 2 2 2 3 2 2 2 3 3 3 3 2 3 3 2 2 2 3 3 2 2 2]
[2 3 2 3 2 2 3 3 2 3 3 3 2 3 2 2 3 3 2 3 2 2 3 2 3 2 2 3 2 2 3 2]
[2 3 2 3 2 3 3 3 3 3 3 2 2 2 2 2 3 2 3 2 2 2 3 2 2 3 3 2 3 2 2 3]
[3 2 3 2 2 3 3 3 2 2 3 3 2 2 3 3 2 3 2 2 3 2 2 3 2 2 3 2 3 2 2 3]
[3 2 3 3 2 3 2 2 2 3 3 3 2 3 3 2 2 2 3 2 3 2 2 3 2 2 3 2 2 2 3 3]]
100 solutions verified out of 133342 x 113193 = 15093381006
我最近一直在算法中使用 numpy 矩阵，但遇到了一个问题：
我一共用了3个矩阵
m1 = [[ 3 2 2 ... 3 2 3]
[ 3 3 3 ... 2 2 2]
[500 501 502 ... 625 626 627]
...
[623 624 625 ... 748 749 750]
[624 625 626 ... 749 750 751]
[625 626 627 ... 750 751 752]]
m1 是 (128,128) 奇异 方阵。前两行看似随机的 2 和 3 序列。接下来的行按照算法从 500 开始填充，从第三行第一列开始为每一行和每一列加一个。
m2 = [[ 2 3 500 ... 623 624 625]
[ 2 2 500 ... 623 624 625]
[ 3 2 500 ... 623 624 625]
...
[ 2 3 500 ... 623 624 625]
[ 2 2 500 ... 623 624 625]
[ 3 2 500 ... 623 624 625]]
m2也是一个(128,128)奇异方阵。这次，随机序列归因于前两列。每行的其余部分填充500、501、502、503等。
m3 = [[ 790 784 157500 ... 196245 196560 196875]
[ 804 811 161000 ... 200606 200928 201250]
[ 180501 180411 36064000 ... 44935744 45007872 45080000]
...
[ 219861 219771 43936000 ... 54744256 54832128 54920000]
[ 220181 220091 44000000 ... 54824000 54912000 55000000]
[ 220501 220411 44064000 ... 54903744 54991872 55080000]]
m3 = m1*m2
所以我想做的是使用 m1 和 m3 恢复 m2。理论上，我所要做的就是执行下面的代码m2 = (m1**-1)*m3。不幸的是，由于 m1 是奇异矩阵，无法计算其逆矩阵，即使可以，矩阵也太大，导致许多数值不精确。
相反，我决定使用 m1 的 Moore-Penrose Inverse，它不需要矩阵是非奇异的，并且与逆矩阵类似，理论上可以恢复 m2, 使用 np.linalg.pinv(m1) * m3.
我再次使用“理论上”一词，因为事实证明，numpy 在涉及大矩阵的此类计算时过于不精确，这是我获得的结果 m2:
[[ 2.46207616 2.48959603 500. ... 623. 624.
625. ]
[ 2.38612549 2.61197086 500. ... 623. 624.
625. ]
[ 2.38711085 2.6125801 500. ... 623. 624.
625. ]
...
[ 2.61998539 2.37184747 500. ... 623. 624.
625. ]
[ 2.54403472 2.4942223 500. ... 623. 624.
625. ]
[ 2.62195611 2.37306595 500. ... 623. 624.
625. ]]
如您所见，m1 的整个 'filler' 部分计算正确，没有问题。但是，前两列似乎有问题，将数字四舍五入为 2s 和 3s 给我一个不正确的 m2.
我正在寻找一种方法，使 np.linalg.pinv() 方法的浮点计算更加精确，以便它可以获得序列的正确值，因为这些值非常重要。
通过一些研究，我了解到 np.linalg.pinv() 有一个名为 rcond 的参数，描述如下：
rcond : (…) array_like of float
Cutoff for small singular values. Singular values smaller (in modulus) than rcond * largest_singular_value (again, in modulus) are set to zero. Broadcasts against the stack of matrices
rcond 默认设置为 1e-15。我认为进一步减少这个数字可能有助于不精确。 1e-16 还不够，从 1e-17 开始，我得到了非常奇怪的值，例如：
[[ 3.000e+00 3.000e+00 5.000e+02 ... 6.230e+02 6.240e+02 6.250e+02]
[ 1.100e+01 4.000e+00 1.840e+02 ... 1.722e+03 2.032e+03 1.831e+03]
[-3.000e+00 -5.000e+00 -4.030e+02 ... -1.232e+03 -7.400e+02 -1.272e+03]
...
[ 1.100e+01 1.200e+01 2.164e+03 ... 4.030e+03 4.872e+03 1.873e+03]
[-1.200e+01 -9.000e+00 -1.618e+03 ... -3.240e+03 -2.519e+03 -4.167e+03]
[ 2.000e+01 2.600e+01 4.535e+03 ... 5.165e+03 5.881e+03 5.189e+03]]
所以，基本上，我被困住了，我不知道如何提高精度。
最糟糕的是，我有一个可以显着提高浮点精度的模块，它称为 mpmath 并且还有一些矩阵似乎更适合我的算法，如 numpy 矩阵。但是 mpmath 没有计算伪逆的方法，并且 numpy 没有将自己的浮点精度调整为 mpmath.
你有什么建议我可以尝试使用伪逆方法得到正确的 m2 吗？
你的烦恼与pinv准确与否无关
正如您自己注意到的那样，您的矩阵严重缺乏秩，m1 的秩为 4 或更低，m2 的秩为 3 或更低。因此，您的系统 m1@x = m3 极度欠定，无法恢复 m2.
即使我们把我们知道的关于 m2 结构的所有信息都扔进去，即前两列 3 和 2，其余 500 向上计数，仍然有大量的组合解。
如果有足够的时间，下面的脚本可以找到它们。在实践中，我没有看到超过 32x32 矩阵，在下面显示的 运行 中产生
15093381006 个不同的有效重建 m2' 满足 m1@m2' = m3 和我刚才提到的结构约束。
import numpy as np
import itertools
def make(n, offset=500):
offset -= 2
res1 = np.add.outer(np.arange(n), np.arange(offset, offset+n))
res1[:2] = np.random.randint(2, 4, (2, n))
res2 = np.add.outer(np.zeros(n, int), np.arange(offset, offset+n))
res2[:, :2] = np.random.randint(2, 4, (n, 2))
return res1, res2
def subsets(L, n, mn, mx, prepend=[]):
if n == 0:
if mx >= mn:
yield prepend
elif n == 1:
for l in L[L.searchsorted(mn):L.searchsorted(mx, 'right')]:
yield prepend + [l]
else:
ps = L.cumsum()
ps[n:] -= ps[:-n]
ps = ps[n-1:]
for i in range(L.searchsorted(mn - np.sum(L[len(L)-n+1:])),
ps.searchsorted(mx, 'right')):
yield from subsets(L[i+1:], n-1, mn - L[i], mx - L[i],
prepend = prepend + [L[i]])
def solve(m1, m3, ci=0, offset=500):
n, n = m1.shape
col = m3.T[ci]
n3s = col[3] - col[2] - 2 * n
six = col[2] - offset * (col[3] - col[2]) - n * (n-1)
idx = np.lexsort(m1[:2])
m1s = m1[:2, idx]
sm = m1s[1].searchsorted(2.5)
sl = m1s[0, :sm].searchsorted(2.5)
sr = sm + m1s[0, sm:].searchsorted(2.5)
n30 = n - sl - sr + sm
n31 = n - sm
n330 = col[0] - 4*n - 2*n3s - 2*n30
n331 = col[1] - 4*n - 2*n3s - 2*n31
for n333 in range(max(0, n330 - sm + sl, n331 - sr + sm),
min(n - sr, n330, n331) + 1):
n332 = n330 - n333
n323 = n331 - n333
n322 = n3s - n332 - n323 - n333
mx333 = six - idx[sl:sl+n332].sum() - idx[sm:sm+n323].sum() \
- idx[:n322].sum()
mn333 = six - idx[sm-n332:sm].sum() - idx[sr-n323:sr].sum() \
- idx[sl-n322:sl].sum()
for L333 in subsets(idx[sr:], n333, mn333, mx333):
mx332 = six - np.sum(L333) - idx[sm:sm+n323].sum() \
- idx[:n322].sum()
mn332 = six - np.sum(L333) - idx[sr-n323:sr].sum() \
- idx[sl-n322:sl].sum()
for L332 in subsets(idx[sl:sm], n332, mn332, mx332,
prepend=L333):
mx323 = six - np.sum(L332) - idx[:n322].sum()
mn323 = six - np.sum(L332) - idx[sl-n322:sl].sum()
for L323 in subsets(idx[sm:sr], n323, mn323, mx323,
prepend=L332):
ex322 = six - np.sum(L323)
yield from subsets(idx[:sl], n322, ex322, ex322,
prepend=L323)
def recon(m1, m3, ci=0, offset=500):
n, n = m1.shape
nsol = nfp = 0
REC = []
for i3s in solve(m1, m3, ci, offset):
rec = np.full(n, 2)
rec[i3s] = 3
if not np.all(m3.T[ci] == m1@rec):
print('!', rec, m3.T[ci], m1@rec)
nfp += 1
else:
nsol += 1
REC.append(rec)
print('col', ci, ':', nsol, 'solutions,', nfp, 'false positives')
return np.array(REC)
def full_recon(m1, m3, offset=500, subsample=10):
n, n = m1.shape
col0, col1 = (recon(m1, m3, i, offset) for i in (0, 1))
yield col0.shape[0], col1.shape[0]
if not subsample is None:
col0, col1 = (col[np.random.choice(col.shape[0], subsample)]
if col.shape[0] > subsample else col
for col in (col0, col1))
print('col 0', col0)
print('col 1', col1)
for c0, c1 in itertools.product(col0, col1):
out = np.add.outer(np.zeros(n, int), np.arange(offset-2, offset+n-2))
out[:, :2] = np.c_[c0, c1]
yield out
def check(m1, m3, offset=500, subsample=10):
for cnt, m2recon in enumerate(full_recon(m1, m3, offset, subsample)):
if cnt == 0:
tot0, tot1 = m2recon
continue
assert np.all(m3 == m1@m2recon)
print(cnt, 'solutions verified out of', tot0, 'x', tot1, '=', tot0 * tot1)
样本运行:
>>> m1, m2 = make(32)
>>> check(m1, m1@m2)
col 0 : 133342 solutions, 0 false positives
col 1 : 113193 solutions, 0 false positives
col 0 [[2 2 3 2 2 2 3 2 3 3 3 3 3 2 3 3 2 2 2 2 2 3 3 3 2 2 2 2 2 2 2 2]
[2 2 3 2 3 3 3 3 3 3 2 2 2 2 2 2 2 2 2 3 3 3 2 2 2 3 2 3 2 2 2 2]
[2 3 3 2 3 3 2 2 2 3 2 3 3 2 2 2 2 3 3 3 2 2 2 2 3 2 2 2 2 2 2 3]
[2 2 3 3 2 2 3 3 3 2 2 3 2 3 3 3 2 2 2 2 2 2 2 3 2 3 3 2 2 2 2 2]
[2 3 3 3 2 3 3 2 2 2 2 3 2 2 3 3 2 2 3 2 2 3 2 2 2 2 2 2 3 3 2 2]
[3 2 2 2 3 2 3 3 3 2 3 2 2 2 2 3 3 3 2 2 2 2 3 3 2 3 2 2 2 2 2 2]
[2 2 2 3 3 3 2 2 3 3 3 3 3 2 2 2 2 3 2 2 2 3 2 2 2 3 2 2 3 2 2 2]
[2 2 2 3 3 2 3 3 3 2 3 2 2 3 3 3 2 3 2 2 2 2 2 2 2 2 2 3 2 3 2 2]
[3 2 3 3 2 2 3 2 2 3 2 3 3 2 2 2 3 2 3 2 2 2 2 3 2 3 2 2 3 2 2 2]
[3 2 2 3 3 2 3 3 2 2 2 3 2 3 2 2 3 2 3 3 2 2 2 2 2 3 2 2 2 2 2 3]]
col 1 [[2 2 2 2 3 3 3 3 3 3 3 3 2 2 3 2 2 2 3 3 2 3 2 2 2 3 2 3 2 2 2 3]
[2 3 3 2 3 3 2 2 2 3 3 3 2 2 3 3 2 3 2 3 2 2 3 2 2 2 2 3 3 2 2 3]
[3 2 3 2 2 3 3 2 2 3 3 3 2 2 3 2 3 2 3 2 3 2 3 3 2 2 2 2 2 3 3 2]
[2 2 3 2 3 3 3 2 3 3 2 3 2 3 2 3 2 3 2 2 3 2 3 2 2 3 2 3 2 2 2 3]
[3 3 2 2 2 2 3 2 3 3 3 3 2 3 3 2 2 3 3 2 2 2 2 2 3 2 2 3 3 3 2 2]
[3 3 2 3 2 2 2 3 3 3 2 2 2 3 2 2 2 3 3 3 3 2 3 3 2 2 2 3 3 2 2 2]
[2 3 2 3 2 2 3 3 2 3 3 3 2 3 2 2 3 3 2 3 2 2 3 2 3 2 2 3 2 2 3 2]
[2 3 2 3 2 3 3 3 3 3 3 2 2 2 2 2 3 2 3 2 2 2 3 2 2 3 3 2 3 2 2 3]
[3 2 3 2 2 3 3 3 2 2 3 3 2 2 3 3 2 3 2 2 3 2 2 3 2 2 3 2 3 2 2 3]
[3 2 3 3 2 3 2 2 2 3 3 3 2 3 3 2 2 2 3 2 3 2 2 3 2 2 3 2 2 2 3 3]]
100 solutions verified out of 133342 x 113193 = 15093381006