1) Похоже, что NumPy не делает различия между вектором-строкой и столбцом.
Там все тоньше. Одномерный массив - это не вектор-строка, а просто одномерный массив, а умножение определено для произвольных многомерных, поэтому позволяет больше вольностей, чем в стандартных матричных операциях. Можно написать
b.reshape(b.shape[0],1) или
b[:,np.newaxis] - и будет вектор-столбец,
b.reshape(1,b.shape[0]) или
b[np.newaxis,:] - вектор-строка. Они будут умножаться только на то, на что положено.
2) Зачем звёздочка перед np.ix_?
np.ix_ поворачивает каждый аргумент вдоль новой оси и возвращает их обратно, то есть его результат в данном случае - пара из вектора-столбца иксов и вектора-строки степеней. А
power нужны два аргумента, вот и приходится распаковывать пару. Можно записать это в виде
x[:,np.newaxis] ** np.arange(4)[np.newaxis,:].
3) Что за значок @ и почему для обозначения всяких матричных умножений не ввели просто * ?
Обычные операции уже много лет действуют поэлементно. Как ни странно, поэлементное умножение нужно довольно часто.