You can use Mata interactively.

A = 1,2\3,4		    // Input a 2 × 2 matrix
b = -1\3		    // Input a 2 × 1 column vector
A'b			    // Compute A'b
A[2,1]			    // Access the (2,1) element of A
A[,2]      		    // Access the second column of A
B = diag(b)		    // Create a 2 × 2 matrix with the elements of b on the diagonal
C = A,B     		    // Combine matrices A and B by column
D = A\B     		    // Or by row
I = I(3)		    // Create a 3 × 3 identity matrix I
trace(A)		    // Compute the trace of A
det(A)			    // And the determinant of A
G*H			    // Multiply G (n × k) by H (k × k)
G:*g  			    // Multiply each column of G (n × k)  by  g (n × 1)
st_data(.,.)                // Load Stata dataset into Mata
st_data((1,5\7,9), "mpg")   // Observations 1 to 5 and 7 to 9 of variable mpg

Or you can write a simple program for regression.

mata
st_view(y=., ., "mpg")
st_view(X=., ., ("weight", "foreign", "cons"))
b = invsym(X'X)*X'y
b
e = y - X*b 
n = rows(X)
k = cols(X)
s2= (e'e)/(n-k)
V = s2*invsym(X'X)
V
se = sqrt(diagonal(V))
(b, se)

Or build a regression system using structures or classes.

class linreg
{
    public:
    
	void		setup()    
	real scalar	N()
        real scalar	k()
        real scalar	K()
        real colvector	b_X()
        real scalar	b_c()
        real colvector	b()
        real matrix	V()
        real scalar	s2()
        real scalar	ee()
        real colvector	yhat()

    protected:
    
        // --------------------------------------- inputs
        pointer(real colvector) scalar  y     
        pointer(real matrix) scalar     X     
        real scalar                     cons  

        // --------------------------------------- derived

        real matrix     XX              // X'X; K x K 
        real colvector  Xy              // X'y; K x K
        real matrix     XXinv           // (X'X)^(-1);   K x K
        real colvector  b               // coefficients; K x 1
        real matrix     V               // variance matrix (VCE); K x K
        real scalar     s2              // variance of e
        real scalar     yy              // y'y
        real scalar     ee              // e'e, error sum of squares
        real scalar     K_adj
        real scalar     k_adj

        /*
            Deviation from mean notation:
                    S = X - mean(X)    X as deviation from mean
                    t = y - mean(y)    y as deviation from mean
        */

        real matrix    SS                // S'S; k x k 
        real colvector St                // S't; K x 1 (sic)
        real matrix    SSinv             // (S'S)^(-1); K x K (sic) 
        real rowvector mean_X            // means of rows of X; 1 x k
        real scalar    mean_y            // mean  of row  of y; 1 x 1
        
        void 		init()
        real matrix	XX()
        real colvector	Xy()
        real scalar	yy()
        real matrix	XXinv()
        real matrix	SSinv()
        void		calc_SSinv_cons()
        real matrix	SS()
        real colvector	St()
        real rowvector	mean_X()
        real scalar	mean_y()
        real scalar	K_adj()
        real scalar	k_adj()
}

void linreg::init()
{
        XX      = .z
        Xy      = .z
        XXinv   = .z 
        b       = .z
        V       = .z
        s2      = .z
        yy      = .z
        ee      = .z

        SS      = .z
        St      = .z
        mean_X  = .z
        mean_y  = .z
        SSinv   = .z
        K_adj   = .z 
        k_adj   = .z 
}


/* -------------------------------------------------- entry points --- */

void linreg::setup( real colvector user_y, 
                    real matrix    user_X,  
                    real scalar    user_cons )
{

        assert(rows(y)==rows(X))

        y    = &user_y
        X    = &user_X
        cons = (user_cons!=0)

        init()              // initialize derived vars

}

real scalar linreg::N() return(rows(*X))
real scalar linreg::k() return(cols(*X))
real scalar linreg::K() return(cols(*X) + cons)


/*
        Throughout this code, distinguish between k and K:

                k = # of independent variables EXCLUDING intercept
                K = # of independent variables INCLUDING intercept

        In models without an intercept, K = k
        In models with    an intercept, K = k + 1

        Note: *X is k x k, X does not include vector of 1s, and 
                 b is K x 1, coefficient include intercept if there is one.

        To code the equivalent of X*b, you extract the first k elements
        of b and add the intercept separately.  To make this easier
        think of

                 b = (b_X, b_c)   

        b_X() returns b_X:  1 x k
        b_c() returns b_c:  1 x 1  containing b[K] or 0. 

        Thus, X*b can be calculated as (*X)*b_X() :+ b_c()
        in all cases.
*/


real colvector linreg::b_X()
{
        if (cons) {         
                if (k())  return( b()[| 1\k() |])
                else      return( J(0, 1, .))
        }
        else              return(b())
}


real scalar linreg::b_c() 
{
        return(cons ? b()[K()] : 0) 
}


real colvector linreg::b()
{
        if (b == .z) {
                if (cons) {
                        b          = SSinv() * St()
                        b[K()]    = mean_y() - mean_X()*b_X()

                }
                else    b = XXinv() * Xy()
        }
        return(b)
}


real matrix linreg::V()
{
        if (V == .z) {
                V = s2() * (cons ? SSinv() : XXinv())
        }
        return(V)
}

real scalar linreg::s2()
{
        if (s2 == .z) {
                s2 = ee() / ( N() - K_adj() )
        }
        return(s2)
}

real scalar linreg::ee()
{
        real colvector   e

        if (ee == .z) {
                e = *y - yhat()
                ee = quadcross(e, e)
        }
        return(ee)
}

real colvector linreg::yhat()
{
        return( (*X)*b_X() :+ b_c() )
}


/* ------------------------------------------internal subroutines --- */

real matrix linreg::XX()
{
        if (XX == .z) {
                XX =  quadcross(*X, cons,  *X, cons)
        }
        return(XX)
}

real colvector linreg::Xy()
{
        if (Xy == .z) {
                Xy =  quadcross(*X, cons,  *y, 0)
        }
        return(Xy)
}

real scalar linreg::yy()
{
        if (yy == .z) {
                yy = quadcross(*y, *y)
        }
        return(yy)
}

real matrix linreg::XXinv()
{
        if (XXinv == .z) {
                XXinv = invsym(XX())
        }
        return(XXinv)
}

real matrix linreg::SSinv()
{
        if (SSinv == .z) {
                if (cons) calc_SSinv_cons()
        }
        return(SSinv)
}


void linreg::calc_SSinv_cons() 
{
        real scalar        K, k, i
        real rowvector        v
        real matrix        kxk, rowK, colK

        K = K()                // # X vars including constant
        k = k()                // # x vars excluding constant

        /* 
            Calculate inverse of S'S = SS().
            This routine called only if cons, meaning K = k + 1.
            SS() is k x k -- it excludes constant. 
            We want the full K x K inverse. 
        */

        // ---------------------------------------- SSinv 1 x 1 case ---
        if (K==1) {
                SSinv = 1/N()
                return
        }

        // ----------------------------------- SSinv 2 x 2 or larger ---

        SSinv = J(K, K, .) // matrix now K x K

                             // Range subscripts:
        kxk  = (1,1 \ k,k)   //   k x k submatrix  rows & cols 1..k
        rowK = (K,1 \ K,k)   //   1 x k submatrix; row K
        colK = (1,K \ k,K)   //   k x 1 submatrix; col K

        SSinv[|kxk |]     = invsym(SS())

        SSinv[|rowK|] = v = -(mean_X() * SSinv[|kxk|])
        SSinv[|colK|] = v'

        // SSinv[ K,K  ]     = 1/N() - mean_X()*SSinv[|colK|]
        SSinv[ K,K  ]     = 1/N() - mean_X()*v'
}


real matrix linreg::SS()
{
        if (SS == .z) {
                if (cons) {
                        SS = quadcrossdev(*X, 0, mean_X(), 
                                          *X, 0, mean_X() )
                }
        }
        return(SS)
}


real colvector linreg::St()
{
        if (St == .z) {
                St = quadcrossdev(*X,      0, mean_X(), 
                                  *y,      0, mean_y()) \ 0
        }
        return(St)
}


real rowvector linreg::mean_X() 
{
        if (mean_X == .z) mean_X = mean(*X)
        return(mean_X)
}


real scalar linreg::mean_y()
{
        if (mean_y == .z) mean_y = mean(*y)
        return(mean_y)
}


real scalar linreg::K_adj()
{
        if (K_adj == .z) {
                if (cons) K_adj = colsum(diagonal(SSinv()) :!= 0)
                else        K_adj = colsum(diagonal(XXinv()) :!= 0)
        }
        return(K_adj)
}


real scalar linreg::k_adj()
{
        if (k_adj == .z) {
                k_adj = k() - (K()-K_adj())
        }
        return(k_adj) 
}

end

Resources

There are many great resources for learning Mata.

The official documentation:

• [M] Mata Reference Manual

We recommend starting with the intro sections:

• Introduction and first session
• Using Mata interactively
• Using Mata with ado-files

The Mata Book: A Book for Serious Programmers and Those Who Want to Be by William Gould

The Mata Matters articles published in the Stata Journal:

• Translating Fortran
• Using views onto the data
• Creating new variables—sounds boring, isn't
• Interactive use
• Precision
• Subscripting
• Structures
• Macros
• File processing
• Overflow, underflow and the IEEE floating-point format

Programming an estimation command in Stata, a series of posts on the Stata Blog.

An Introduction to Mata from the Social Science Computing Cooperative at the University of Madison–Wisconsin

Estimation of the Mixture of Gaussians Using Mata, blog post by Dallakyan Aramayis.

Learn more about all of Mata's features.