what do we actually study in machine learning?
I want to point out one aspect which is not adequately reflected in the previous answers. When you start with linear regression, you basically I) propose some model (a way how to represent the relation between $y$ and $X$), e.g. a linear function, and II) you specify a loss function (the sum of squared residuals). The aim is to find coefficients (or weights), which map $X$ to $y$ in a way so that there is no “better way” to predict $y$ (viz. minimise the sum of square residuals).

In the case of linear regression, you have two options: solve the problem analytically or by gradient decent. Optimal coefficients/weights (vector beta) are given by:

$$ (X'X)^{-1} X'y=\beta. $$

So there is no need for gradient decent in this case. It is just simple matrix algebra. 

Here is a minimal example in R, demonstrating how it works:

Source: https://www.r-bloggers.com/regression-via-gradient-descent-in-r/

# Input some data
x0 &lt;- c(1,1,1,1,1) 
x1 &lt;- c(1,2,3,4,5)
x2 &lt;- c(8,4,3,1,8)
x &lt;- as.matrix(cbind(x0,x1,x2))
y &lt;- as.matrix(c(3,7,5,11,14))

# (X'X)^-1 X'y
beta1 = solve(t(x)%*%x) %*% t(x)%*%y 

# R's regression command
beta2 = summary(lm(y ~ x[, 2:3]))

# Gradient decent
m &lt;- nrow(y)
grad &lt;- function(x, y, theta) {
  gradient &lt;- (1/m)* (t(x) %*% ((x %*% t(theta)) - y))
  return(t(gradient))
}

# Define gradient descent update algorithm
grad.descent &lt;- function(x, maxit){
  theta &lt;- matrix(c(0, 0, 0), nrow=1) # Initialize the parameters

  alpha = .05 # Learning rate
  for (i in 1:maxit) {
    theta &lt;- theta - alpha  * grad(x, y, theta)   
  }
  return(theta)
}

# Results
print(grad.descent(x,2000))
print(beta1)
print(beta2)


Results:

&gt; print(grad.descent(x,2000))
            x0     x1        x2
[1,] -0.672238 2.6459 0.1530299
&gt; print(beta1)
         [,1]
x0 -0.6722955
x1  2.6459103
x2  0.1530343
&gt; print(beta2)

Call:
lm(formula = y ~ x[, 2:3])

Residuals:
      1       2       3       4       5 
-0.1979  1.7683 -2.7245  0.9356  0.2185 

Coefficients:
            Estimate Std. Error t value Pr(&gt;|t|)  
(Intercept)  -0.6723     3.3565  -0.200   0.8598  
x[, 2:3]x1    2.6459     0.7677   3.447   0.0748 .
x[, 2:3]x2    0.1530     0.3897   0.393   0.7325  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.399 on 2 degrees of freedom
Multiple R-squared:  0.8561,    Adjusted R-squared:  0.7122 
F-statistic: 5.949 on 2 and 2 DF,  p-value: 0.1439


So as you can see, there is not much of a machine involved in solving the problem (apart of doing the math of course). So why machine learning? Increasing computational power has made it possible to solve very complex problems using gradient decent. The idea is that you depart from simple solutions and look for better, more complex ways to map $X$ to $y$.

You can think about changing the base procedure, e.g. use decision trees instead of a linear function. Since one tree may be only a weak learner, you can combine (or ensemble) trees to a random forest. Or you go even further and use boosting, in which you repeatedly use information contained in the residual to improve your fit. Usually very small trees are used to do this, but you can boost on linear regression as well. 

Here is a little example for L2 boosting:

(from: Boosting Algorithms: Regularization, Prediction and Model Fitting)

library(glmnet)

# Simple boosting example from:
# P. Bühlmann, T. Hothorn (2007), "Boosting Algorithms: Regularization, Prediction and Model Fitting", Statistical Science 22(4), p. 477-505.
# L2-Boosting from Section 3.3 (p. 483)

#######################
# Data 
x = matrix(c(130,165,150,150,140,198,307,350,318,304,302,429),nrow=6,ncol=2,byrow = TRUE) 
y = c(18,15,18,16,17,15)

# Parameter
# This works like a learning rate
nu = 0.1
# Early stopping (SSR value to stop)
es = 7.58765583743327
# Max. iterations (of boosting)
maxiter = 10000

#######################
# OLS benchmark model
olsdata = data.frame(cbind(y,x))
ols = lm(y~.,data=olsdata)
olspreds = predict(ols,newdata=olsdata)
olsssr = sum((y-olspreds)^2)
print(paste0("OLS SSE: ",olsssr))

#######################
# Boosting 

# Initialize f0
f0 = mean(y)
f0ssr = sum((y-f0)^2)

# Lists to store results per boosting iteration
ssrlist = list()
bstep = list()

# Boosting (p. 483, Sec. 3.3, L2-Boosting)
for (n in seq(1:maxiter)){
  # I) Get residual
  if(n==1){f=f0} # initialize in first step  
  u=y-f          # get residual from estimation 
  # II) Fit residual / predict, I use ridge (alpha=0)
  reg = glmnet(x,u,alpha=0, lambda=2)       
  g = predict(reg, newx=x, type="link") 
  # III) Update f
  f=f+nu*g
  # Print feedback
  cat(paste0("Step: ", n," SSR: ", sum(u^2), "\n" ))
  # Save SSR/iter to list
  ssrlist[[n]] &lt;- sum(u^2)
  bstep[[n]] &lt;- n
  # Early stopping rule
  if(sum(u^2)&lt;es){break}
}

#######################
# End statement
cat("\n ~~~ RESULTS ~~~ \n")
cat(paste0("OLS SSR:           ", olsssr, "\n"))
cat(paste0("Initialised SSR:   ", f0ssr, "\n"))
cat(paste0("Last SSR boosting: ", ssrlist[length(ssrlist)], " [after ",bstep[length(bstep)], " iterations] \n"))


The problem here is that you do not only need to solve one regression, but 10000 (in the example here). In real world applications, the amount of math that has to be done is immense (counts for boosting and neural nets as well).

So what is machine learning? You want to let “the machine learn things” (recognise patterns in data, recognise images, make prediction based on that etc) and you will likely be unable to do this learning task on your own using pen and paper. Improved computing capacity over time has enabled us to solve very complex problems. And because of this cheap resource, new tools and methods develop. You can do things today on your laptop which literally seemed impossible in the late 1990s.

Of course you need to tell “the machine” how to learn. This is a lot of statistics but also computer science. Neural nets for example have mostly been advanced by computer scientists. Ultimately, the machine is merely a tool, but a very useful one.