In the FEM we usually use a trick to find the coefficients for shape functions:
Finding $M^{-1}$ in the system $MC=I$, where for example in the linear case: $$ \begin{aligned} M&=[1 ,x_1,y_1;1 ,x_2,y_2;1 ,x_3,y_3],\\ C&=[a_1, b_1,c_1;a_2, b_2,c_2;a_3, b_3,c_3],\\ I&=[1,0,0;0,1,0;0,0,1]. \end{aligned} $$ and $\phi_1=a1+b_1x_1+c_1y_1$ and also $\phi_2,\phi_3$ are defined in the same way. Since $MC=I$ so $M$ is inverse of $C$. But how do we ensure that $M$ is invertible for a higher degree basis?
This is ensured by choosing the degrees of freedom appropriately. (That's why you need, e.g., additional degrees of freedom on the edges of a triangle for quadratic elements.)
Basically, constructing a (local) finite element space consists of choosing a reference element $K$ (e.g., a triangle) and a polynomial space $P$ (e.g., bivariate polynomials on $K$ of total order up to $2$). You then need to find a set of conditions (e.g., point evaluation in some $x_i\in K$) to make the interpolation problem "for given $f\in C(K)$, find $p\in P$ such that $p(x_i) = f(x_i)$ for all $i$" have a unique solution -- these are then the degrees of freedom you need on each element.
It is then a theorem of linear algebra that the matrix $M$ (which is called Vandermonde matrix) is always invertible. (Note that your equation $MC=I$ is just a compact way of writing the nodal interpolation conditions $\phi_i(x_j,y_j) = \delta_{ij}$.)
