matmul.cpp потвърждава, че PCA::Operator()
се използва от PCACompute()
, но собствените стойности се отхвърлят. Така че направих това:
# The following mimics PCA::operator() implementation from OpenCV's
# matmul.cpp() which is wrapped by Python cv2.PCACompute(). We can't
# use PCACompute() though as it discards the eigenvalues.
# Scrambled is faster for nVariables >> nObservations. Bitmask is 0 and
# therefore default / redundant, but included to abide by online docs.
covar, mean = cv2.calcCovarMatrix(PCAInput, cv2.cv.CV_COVAR_SCALE |
cv2.cv.CV_COVAR_ROWS |
cv2.cv.CV_COVAR_SCRAMBLED)
eVal, eVec = cv2.eigen(covar, computeEigenvectors=True)[1:]
# Conversion + normalisation required due to 'scrambled' mode
eVec = cv2.gemm(eVec, PCAInput - mean, 1, None, 0)
# apply_along_axis() slices 1D rows, but normalize() returns 4x1 vectors
eVec = numpy.apply_along_axis(lambda n: cv2.normalize(n).flat, 1, eVec)
(опростяване на предположенията: редове = наблюдения, колони = променливи; и има много повече променливи от наблюдения. И двете са верни в моя случай.)
Това до голяма степен работи. По-долу old_eVec
е резултатът от cv2.PCACompute()
:
In [101]: eVec
Out[101]:
array([[ 3.69396088e-05, 1.66745325e-05, 4.97117583e-05, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -7.23531536e-06, -3.07411122e-06, -9.58259793e-06, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 1.01496237e-05, 4.60048715e-06, 1.33919606e-05, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
...,
[ -1.42024751e-04, 5.21386198e-05, 3.59923394e-04, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -5.28685812e-05, 8.50139472e-05, -3.13278542e-04, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 2.96546917e-04, 1.23437674e-04, 4.98598461e-04, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])
In [102]: old_eVec
Out[102]:
array([[ 3.69395821e-05, 1.66745194e-05, 4.97117981e-05, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -7.23533140e-06, -3.07411415e-06, -9.58260534e-06, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 1.01496662e-05, 4.60050160e-06, 1.33920075e-05, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
...,
[ -1.42029530e-04, 5.21366564e-05, 3.60067672e-04, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -5.29163444e-05, 8.50261567e-05, -3.13150231e-04, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -7.13724992e-04, -8.52700090e-04, 1.57953508e-03, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]], dtype=float32)
Има някаква загуба на прецизност, видима към края на резултатите (въпреки че в действителност бързото начертаване на абсолютната разлика не разкрива модел на неточността).
57% от елементите имат различна от нула абсолютна разлика.
От тях 95% се различават с по-малко от 2e-16 и средната AD е 5,3e-4 - въпреки това AD може да достигне до 0,059, което е много, когато вземете предвид, че всички стойности на собствения вектор са между -0,048 до 0,045.
Има код в PCA::Operator()
, който преобразува в най-големия ctype; от друга страна old_eVec
беше float32 в сравнение с моя собствен код, произвеждащ float64. Струва си да се спомене, че при компилирането на numpy получих някои грешки, свързани с точността.
Като цяло загубата на прецизност изглежда е свързана с ниски собствени вектори със собствени стойности, което отново сочи към грешка при закръгляване и т.н. Горната реализация дава резултати, подобни на PCACompute(), като дублира поведението.
person
benxyzzy
schedule
24.05.2015