Требуется немного творческой индексации массивов, но, надеюсь, это то, что вам нужно. Примечание: в данном случае R не является диагональной матрицей.
import numpy
A = numpy.array([[13865.995, 14020.785, 14020.788 ],
[ 6885.634, 6784.8813, 6837.428 ],
[14510.425, 14712.554, 14657.217 ],
[ 7688.0923, 7817.8457, 7792.413 ],
[10473.903, 10417.55, 10469.508 ],
[10485.661, 10348.152, 10632.414 ]])
B = numpy.array([[14662.705, 14869.951, 15166.294 ],
[13865.995, 14020.785, 14020.788 ],
[ 9780.559, 10038.395, 10202.31 ],
[ 6885.634, 6784.8813, 6837.428 ],
[ 7167.9575, 7357.9062, 7287.3003],
[14510.425, 14712.554, 14657.217 ],
[12825.017, 12680.751, 12823.563 ],
[ 7688.0923, 7817.8457, 7792.413 ],
[ 6861.9443, 6826.6245, 6758.8965],
[10473.903, 10417.55, 10469.508 ],
[ 8498.976, 8637.245, 8718.052 ],
[10485.661, 10348.152, 10632.414 ]])
a = numpy.reshape(A,(A.size,1))
b = numpy.reshape(B,(B.size,1))
R = numpy.zeros((A.size,B.size))
R[[a for a in range(A.size)],[a for a in range(B.size) if (a // 3) % 2]] = 1
print( numpy.all(numpy.dot(R, b) == a) )
# True