For most appliances, both will give you the same results. Generally the standard pythonic a*a or a**2 is faster than the numpy.square() or numpy.pow(), but the numpy functions are often more flexible and precise. If you do calculations that need to be very accurate, stick to numpy and probably even use other datatypes float96.

For normal usage a**2 will do a good job and way faster job than numpy. The guys in this thread gave some good examples to a similar questions.

>>> import numpy
>>> print numpy.power.__doc__

power(x1, x2[, out])

First array elements raised to powers from second array, element-wise.

Raise each base in `x1` to the positionally-corresponding power in
`x2`.  `x1` and `x2` must be broadcastable to the same shape.

Parameters
----------
x1 : array_like
    The bases.
x2 : array_like
    The exponents.

Returns
-------
y : ndarray
    The bases in `x1` raised to the exponents in `x2`.

Examples
--------
Cube each element in a list.

>>> x1 = range(6)
>>> x1
[0, 1, 2, 3, 4, 5]
>>> np.power(x1, 3)
array([  0,   1,   8,  27,  64, 125])

Raise the bases to different exponents.

>>> x2 = [1.0, 2.0, 3.0, 3.0, 2.0, 1.0]
>>> np.power(x1, x2)
array([  0.,   1.,   8.,  27.,  16.,   5.])

The effect of broadcasting.

>>> x2 = np.array([[1, 2, 3, 3, 2, 1], [1, 2, 3, 3, 2, 1]])
>>> x2
array([[1, 2, 3, 3, 2, 1],
       [1, 2, 3, 3, 2, 1]])
>>> np.power(x1, x2)
array([[ 0,  1,  8, 27, 16,  5],
       [ 0,  1,  8, 27, 16,  5]])
>>>

Precision

As per the discussed observation on numerical precision as per @GarethRees objection in comments:

>>> a = numpy.ones( (3,3), dtype = numpy.float96 ) # yields exact output
>>> a[0,0] = 0.46002700024131926
>>> a
array([[ 0.460027,  1.0,  1.0],
       [ 1.0,  1.0,  1.0],
       [ 1.0,  1.0,  1.0]], dtype=float96)
>>> b = numpy.power( a, 2 )
>>> b
array([[ 0.21162484,  1.0,  1.0],
       [ 1.0,  1.0,  1.0],
       [ 1.0,  1.0,  1.0]], dtype=float96)

>>> a.dtype
dtype('float96')
>>> a[0,0]
0.46002700024131926
>>> b[0,0]
0.21162484095102677

>>> print b[0,0]
0.211624840951
>>> print a[0,0]
0.460027000241

Performance

>>> c    = numpy.random.random( ( 1000, 1000 ) ).astype( numpy.float96 )

>>> import zmq
>>> aClk = zmq.Stopwatch()

>>> aClk.start(), c**2, aClk.stop()
(None, array([[ ...]], dtype=float96), 5663L)                #   5 663 [usec]

>>> aClk.start(), c*c, aClk.stop()
(None, array([[ ...]], dtype=float96), 6395L)                #   6 395 [usec]

>>> aClk.start(), c[:,:]*c[:,:], aClk.stop()
(None, array([[ ...]], dtype=float96), 6930L)                #   6 930 [usec]

>>> aClk.start(), c[:,:]**2, aClk.stop()
(None, array([[ ...]], dtype=float96), 6285L)                #   6 285 [usec]

>>> aClk.start(), numpy.power( c, 2 ), aClk.stop()
(None, array([[ ... ]], dtype=float96), 384515L)             # 384 515 [usec]