In my previous blog, the focus was on building an Ordinary Least Squares (OLS) regression model using SAS's algebraic modeling language, the OPTMODEL procedure, which is designed for building and solving optimization models. In this blog, we'll focus on another type of linear regression model called Least Absolute Deviation (LAD).   Unlike OLS regression, which minimizes the sum of the squared residuals, LAD regression minimizes the sum of the absolute value of the residuals. LAD is considered a more robust regression model since residuals are penalized linearly (as opposed to quadratically in OLS).   The tradeoff for increased robustness is a larger residual variance, assuming the error distribution is normal. Depending on the circumstances, LAD is often a preferred alternative to OLS when known outliers are present in the data.   Within the model formulation, the only difference between LAD and OLS is in the objective function.   The objective function for LAD is shown below, which minimizes the sum of the absolute value of the residuals:   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   In a perfect world, we could simply modify the objective function from the OLS formulation, re-run the model, and output the new results. However, the problem is that LAD regression is more computationally difficult to solve than OLS due to the presence of the absolute value function, which introduces non-smoothness into the optimization problem. This "non-smoothness", for lack of a better term, makes it difficult to find a closed-form solution, often resulting in convergence issues for the optimization algorithm, causing it to iterate indefinitely and ultimately failing to find the optimal solution.   But don't take my word for it. Let's try it out.   Below are the set, parameter, and decision variable declarations copied over from the OLS formulation in my previous blog, along with the read data statement, which reads in the sashelp.baseball data set into the OPTMODEL procedure and populates the sets and parameters accordingly. proc optmodel; set <str> PLAYERS; set IVARS = /'nHits' 'nBB'/; num y{PLAYERS}; num x{PLAYERS,IVARS}; read data sashelp.baseball into PLAYERS=[Name] y=nRuns {k in IVARS} <x[Name,k]=col(k)>; var Intercept; var Beta{IVARS};   The new LAD objective function below is the only modification made from the original OLS model in the previous blog. SAS programmers will recognize the abs() function is being used to represent absolute value.   min Obj = sum{i in PLAYERS} abs(y[i] - (sum{k in IVARS} Beta[k]*x[i,k] + Intercept));   The solve statement is specified below. The create data statement creates an output data set of the observed and predicted values for each player. Recall that when you specify the default or standalone solve statement, OPTMODEL will determine and apply the appropriate solver from a set of optimization solvers including LP, MILP, NLP, and QP.   solve; create data work.lad_output from [player] = {i in PLAYERS} y pred=(Intercept + sum{k in IVARS} Beta[k]*x[i,k]); quit;   The results below show the NLP solver was called, and due to the reasons mentioned above, the iteration limit was reached without arriving at an optimal solution.     Well that didn't work! What do we do now? Well... we increase the iteration limit of course! From the Solution Summary table, the default maximum iteration limit is shown to be 5,000. Let's boost that up to 100,000! To do this, I'll modify the solve statement by calling the NLP solver explicitly and setting the maxiter= option to 100,000.   solve with nlp/maxiter=100000;     Both the log and Solution Summary table indicate the model failed to converge. Interestingly, the process suspended after 53 iterations, suggesting that simply increasing the number of iterations isn't going to help solve this problem.   At this point, you may be tempted to throw up your hands, say "to heck with it", and just use OLS instead. But stay with me. There's an easy workaround.     In fact, there's another solver option within the OPTMODEL procedure that will quickly and easily solve LAD to optimality. That option is the solve linearize option.    Linearization in optimization is a powerful technique used to reformulate nonlinear functions with simpler linear functions, allowing us to solve previously unsolvable problems. By replacing complex nonlinear relationships (e.g., absolute value, etc.) with its linear equivalent formulation, optimization algorithms can be efficiently applied to solve the problem.   In other words, linearization simplifies the optimization process by transforming nonlinear objective functions (or constraints) into their linear equivalents, making them amenable to the analytical and numerical techniques required for finding optimal solutions.   The solve linearize statement in the OPTMODEL procedure contains numerous linearization techniques to automatically reformulate nonlinear models into their linear equivalent form. This means that we, as optimization practitioners, do not need to be experts or memorize all of the linearization techniques to easily take advantage of them!    Performing automatic linearization of our LAD model is as simple as modifying the solve statement and re-running the program:   solve linearize;   The Solution Summary table now shows the linear programming (LP) solver was called using the default Dual Simplex algorithm. For the first time, we see the solution status is Optimal with a solution time of 0.02 seconds.     Printing the decision variables (i.e., parameter estimate values) to the Results Viewer shows the new model estimates.     Below is a comparison to the OLS estimates from the previous blog:     OLS LAD Intercept -4.233 -4.861 nBB 0.3008 0.3077 nHits 0.4299 0.4336   Problem solved. The solve linearize statement automatically reformulated the nonlinear LAD model into its linear equivalent and quickly solved it to optimality. The output data set provides us with the observed and predicted values for each player in the sashelp.baseball data set, allowing us to further construct any of the desired plots and model diagnostics of our choosing.   As a result, I had originally intended to conclude the blog here. However, I'd argue that half the fun is learning what the solve linearize statement did behind the scenes that allowed our LAD model to solve to optimality.   Quick side story: I used to love visiting magic shops as a kid. I was always enamored with all of the seemingly impossible magic tricks performed by the owner to entice us to purchase the trick. They would perform the magic trick for you, but if you wanted to learn how it worked, you had to purchase it. Only after you purchased it would they show you the secret(s) and teach you to perform it yourself. Even to this day, I love visiting magic shops, and The Prestige is, still to this day, one of my all time favorite movies.     However, unlike magic tricks, the solve linearize statement has nothing to do with slight of hand, hidden compartments, or smoke and mirrors. In fact, if we knew the linearization technique off the top of our heads, we could program the reformulated model directly inside of the OPTMODEL procedure and bypass the solve linearize statement altogether. But for now let's assume we don't.    In order to reverse engineer what this one instance of solve linearize did for us automatically, we first need to determine the new model formulation that was used. To do this, you can specify expand/linearize; directly after the solve linearize; statement: solve linearize; expand / linearize;   This will write the reformulated linear model to the Results Viewer.   ** Pro Tip **: When doing this, first modify the size of your problem.   To keep the expand/linearize; output manageable, the read data statement is modified to include only players from the Atlanta Braves with last names beginning with H or M. This effectively reduces the set of players from 322 to 5. The put statement writes each of the five elements in the PLAYERS set to the log.   read data sashelp.baseball(where=(Team='Atlanta' and substr(Name,1,1) in ('H','M'))) into PLAYERS=[Name] y=nRuns {k in IVARS} <x[Name,k]=col(k)>; put PLAYERS=;       Running the solve linearize statement now with the expand/linearize option on the reduced set of players shows the linear reformulation that the OPTMODEL procedure performed behind the scenes.     The expand/linearize output shows the OPTMODEL procedure added a set of decision variables for each player: _ADDED_VAR_[1], [2], ... , [5]. These are the absolute values of the residual for each player.   Two sets of constraints are added to the model for each player, _ADDED_CON_[1], [2], ... , [10], and serve as bounds to constrain the added set of decision variables depending on if the observed value is greater than or less than the predicted value.   Let's now replicate this new model formulation using the OPTMODEL procedure.   First I'll create a new set of decision variables, Z, for each player, corresponding to _ADDED_VAR_[1], [2], ... , [5].    var Z{PLAYERS}; According to the expand/linearize; output, the new objective function minimizes the sum of Z across all elements in the PLAYERS set.   min Obj = sum{i in PLAYERS} Z[i]; For the two sets of constraints in the expand/linearize; output, notice that for each player that the linear combination of decision variables all reside on the left-hand side of the inequality. This is done intentionally, but from an applied modeling perspective, we can re-write these constraints any way we'd like to, so long as they're algebraically equivalent to the ones above.   con y_larger{i in PLAYERS}: Z[i] >= y[i] - (sum{k in IVARS} Beta[k]*x[i,k] + Intercept); con yhat_larger{i in PLAYERS}: Z[i] >= - y[i] + (sum{k in IVARS} Beta[k]*x[i,k] + Intercept);   Take the first player in the reduced set, Bob Horner, as an example. Bob's observed value (nRuns) is 70, denoted as y[i] in the OPTMODEL code. Suppose the model predicts a value for Bob less than 70. Let's say it predicts 65. The residual would be +5 (i.e., 70 - 65).   In the first constraint, Z[i] >= 5, and in the second constraint, Z[i] >= -5. The smallest value Z[i] is allowed to take while satisfying both constraints is 5, which is the absolute value of Bob's residual.   If the model instead predicts a value for Bob greater than 70, say 80, then his residual is -10 (i.e., 70 - 80). In the first constraint, Z[i] >= -10, and in the second constraint, Z[i] >= 10. Again, the smallest value Z[i] is allowed to take while satisfying both constraints is 10, which is the absolute value of Bob's residual.   This new model formulation is the linear equivalent to the original formulation containing the nonlinear objective function, and with the help of the expand/linearize; statement, we've effectively reverse engineered what the solve linearize statement did automatically behind the scenes for us.   I hope you found this blog useful and intuitive, and I encourage you to explore all of the automatic linearization techniques packed into the solve linearize; statement. Others include minimax and maximin objective functions, the product of two binary variables, the product of a binary and bounded variable, etc.    See below for the complete list of programs discussed above and try them out for yourself!   --------------------------------------------------------------------------------------------------------------------------   Nonlinear (original) model formulation with and without the solve linearize; and expand/linearize; statements: proc optmodel; set <str> PLAYERS; set IVARS = /'nHits' 'nBB'/; num y{PLAYERS}; num x{PLAYERS,IVARS}; read data sashelp.baseball into PLAYERS=[Name] y=nRuns {k in IVARS} <x[Name,k]=col(k)>; *read data sashelp.baseball(where=(Team='Atlanta' and substr(Name,1,1) in ('H','M'))) into PLAYERS=[Name] y=nRuns {k in IVARS} <x[Name,k]=col(k)>; var Intercept; var Beta{IVARS}; min Obj = sum{i in PLAYERS} abs(y[i] - (sum{k in IVARS} Beta[k]*x[i,k] + Intercept)); solve; *solve with nlp/maxiter=100000; *solve linearize; *expand / linearize; create data work.lad_output from [player] = {i in PLAYERS} y pred=(Intercept + sum{k in IVARS} Beta[k]*x[i,k]); quit;   Linear formulation:   proc optmodel; set <str> PLAYERS; set IVARS = /'nHits' 'nBB'/; num y{PLAYERS}; num x{PLAYERS,IVARS}; read data sashelp.baseball into PLAYERS=[Name] y=nRuns {k in IVARS} <x[Name,k]=col(k)>; var Intercept; var Beta{IVARS}; var Z{PLAYERS}; min Obj = sum{i in PLAYERS} Z[i]; con y_larger{i in PLAYERS}: Z[i] >= y[i] - (sum{k in IVARS} Beta[k]*x[i,k] + Intercept); con yhat_larger{i in PLAYERS}: Z[i] >= - y[i] + (sum{k in IVARS} Beta[k]*x[i,k] + Intercept); solve; create data work.output from [player] = {i in PLAYERS} y z pred=(Intercept + sum{k in IVARS} Beta[k]*x[i,k]); quit; ... View more

  This purpose of this series is to demystify the optimization model building process by building several familiar linear regression models in SAS using the OPTMODEL procedure. This post is geared towards data scientists familiar with linear regression models who are interested in learning more about mathematical programming.   The easiest way to learn a new data science topic (e.g., mathematical programming), in my opinion, is to anchor it to an existing topic that you're already familiar with, so this post assumes you understand the fundamentals of linear regression.   The data being used are the sashelp.baseball data set, which contains 322 players. The dependent variable (y) is nRuns, which is the number of runs a player scored during the season. The two independent variables (x1, x2) being used are nHits and nBB, which are the number of hits and walks, respectively, for each player during the season.   The screen shot below shows the raw data for the first 15 players.     Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   There are numerous procedures and techniques for running ordinary least squares regression in SAS. If you prefer to do it in matrix algebra, the IML procedure is an excellent resource. Other procedures in SAS for OLS are REG, GLMSELECT, GLM, etc. The REG procedure code and output is below.   proc reg data=sashelp.baseball; model nRuns = nHits nBB; run; quit;   The general functional form of the model is provided below.     Even with only two independent variables, the linear predictor appears to do an adequate job fitting the data.       From the Parameter Estimate column above, the linear predictor is constructed from the general functional form:     Recall that in OLS, the objective is to find parameter estimate values that minimize the residual sum of squares.     Or written the "long" way:     Replacing the actual parameter estimates from the model with the Betas in the equation above, the output from the REG procedure effectively told us this function is minimized with the following parameter estimate values:     If you're still tracking with me to this point, we have everything we need to reconstruct the OLS model as a quadratic programming problem within the OPTMODEL procedure.   A mathematical programming problem can be broken down into basic components that are used to solve the problem. It's similar to baking a cake, where the cake can be broken down into the various ingredients used to make it. Similarly, a mathematical programming problem can be broken into its own "ingredients" used to formulate it.   The first "ingredient" are sets. A set is a collection of distinct objects (or elements) that share common characteristics. For this problem, we have two sets: a set of observations (i.e., players), and a set of independent variables (i.e., nHits and nBB).   Inside of the OPTMODEL procedure, sets can be read from SAS data sets or defined explicitly from within the procedure. For example, the code below declares a character set of PLAYERS. When you hear the word "declare", think "create a memory space for". In OPTMODEL, when the data come from external data sources, you must declare sets before populating them with values.   proc optmodel; set <str> PLAYERS;   The Name column from sashelp.baseball above contains one distinct player per row. In mathematical programming, each distinct player name is an element in the set named PLAYERS. The set name is user-defined, and "PLAYERS" seems to be a reasonable name for the elements it will contain. If you recall, there are 322 players in the sashelp.baseball data set, meaning there will be 322 elements in the PLAYERS set once the data has been read in. The "str" within angle brackets indicates that PLAYERS is a character set, since the Name column is character. If it were instead a numeric set, you would instead specify "num" within angle brackets.   The second set is a space-delimited list of independent variables called IVARS. Unlike the PLAYERS set, each element in IVARS is explicitly defined between a pair of forward slashes "/". When explicitly defining elements within a set, "str" or "num" within angle brackets is not required, as the OPTMODEL procedure is able to infer the type from the hard-keyed elements.   set IVARS = /'nHits' 'nBB'/;   The next ingredient are parameters. Do not confuse parameters with parameter estimates! A parameter in a mathematical programming model is simply a known constant that is used throughout the model. In other words, what input data are we using to build the model? We have our known target variable nRuns for each player, along with our known independent variables for each player, nHits and nBB.   These are the three numeric parameters we'll create inside of the OPTMODEL procedure. To do this, we'll use the num statement.   The first parameter, y, corresponds to the target variable nRuns, and is indexed by PLAYERS using curly brackets { }. Think of it like an Nx1 vector, where y[1] is the number of runs scored by the first player in the PLAYERS set, y[2] is the number of runs scored by the second player in the PLAYERS set, and so on.   num y{PLAYERS};   The next parameter, x, is indexed by both the PLAYERS set and the IVARS set. Think of it like a 322x2 matrix with players down the rows and the independent variables along the columns.   num x{PLAYERS,IVARS};   The read data statement reads the sashelp.baseball data set into the OPTMODEL procedure and populates the PLAYERS set and parameters accordingly.   read data sashelp.baseball into PLAYERS=[Name] y=nRuns {k in IVARS} <x[Name,k]=col(k)>;   For a sanity check, you can now optionally write the elements in the PLAYERS set to the log by typing:   put PLAYERS=;   Additionally, you have the option to print the parameters to the Results Viewer.   print y; print x;   The third ingredient in a mathematical programming model are the decision variables. The decision variables are the unknown variables that you want to find the optimal values of. In our OLS regression model, these are the betas.   In mathematical programming, decision variables are allowed to take on one of three types: continuous, integer, or binary. For continuous and integer decision variables, user-defined lower and/or upper bounds can be applied. For OLS, the betas are unbounded and allowed to take fractional values.   In this example, we have three decision variables: B0 (intercept), B1, and B2.   Decision variables are declared using the var statement in OPTMODEL. The first decision variable, Intercept, corresponds to B0 in the functional form above. The set of decision variables Beta are indexed by the IVARS set, meaning the OPTMODEL procedure will create as many Beta decision variables as there are elements in the IVARS set. Since there are only two elements in IVARS, "Beta{IVARS}" will create one decision variable for nHits and one for nBB.   var Intercept; var Beta{IVARS};   The fourth ingredient in a mathematical programming model is the objective function. In the world of machine learning, it's sometimes referred to as the loss function, implying a minimization problem, but the objective function in a mathematical programming model can be either a maximization or a minimization function.   Using the sum{ } aggregation operator to represent big Sigma, the objective function in the OPTMODEL procedure mimics the functional form below, minimizing the residual sum of squares.     min Obj = sum{i in PLAYERS} (y[i] - (Intercept + sum{k in IVARS} Beta[k]*x[i,k]))**2;   To turn it into a maximization problem, you could simply multiply the objective function above by -1.   max Obj = -1*(sum{i in PLAYERS} (y[i] - (Intercept + sum{k in IVARS} Beta[k]*x[i,k]))**2);   The last ingredient in a mathematical programming model are constraints. Constraints are specific conditions or limitations that restrict the feasible solutions in an optimization problem. There are no constraints in the OLS model. In mathematical programming, this is called an unconstrained model.   The term "linear" regression means the model is linear with respect to the Betas, or decision variables. In statistics you may've heard it as "linear in the parameters", which is also referring to the Betas, but as mentioned above, the term "parameter" can be confusing since it means different things to different audiences.   To solve this using mathematical programming, we need to use a quadratic programming solver (not a linear programming solver), due to the squared term in the objective function.   For most optimization models, you can simply specify "solve;" and the OPTMODEL procedure will automatically determine and apply the appropriate solver for your model.   solve;   The remaining code below prints the optimal decision variable values to the Results Viewer, and outputs a SAS data set called work.output containing each player, along with the observed and predicted value for the target variable nRuns.   print Intercept; print Beta; create data work.output from [player] = {i in PLAYERS} y pred=(Intercept + sum{k in IVARS} Beta[k]*x[i,k]); quit;   Among other things, the Solution Summary table confirms the QP solver was used, specifically the Interior Point algorithm. The Solution Status indicates the solution is optimal, the objective function is minimized at 25,853.92, and the solution time was 0.01 seconds.     The optimal decision variable values (i.e., the Betas) match the output from the REG procedure above.     Lastly, the first 15 observations from the work.output table are shown below, including both the observed and predicted values for the target variable nRuns.     In a future post, we'll use the OPTMODEL procedure to solve another linear programming model called Least Absolute Deviation, or LAD regression. It's similar to OLS, but instead of minimizing the residual sum of squares, the objective function in LAD minimizes the sum of the absolute deviations between the observed and predicted values. LAD is considered a more robust regression model since residuals are penalized linearly (as opposed to quadratically in OLS).   In summary, the purpose of this post is to show how OLS can be constructed as a mathematical programming model. The OPTMODEL procedure is classified as an algebraic modeling language which allows us to build and solve optimization models of varying size and complexity. For those interested in learning more about how to build mathematical programming models using the OPTMODEL procedure but are wondering where to start, hopefully this example provided you with a gentle introduction to the process along with a straightforward introduction to the programming syntax. For another introductory example, check out how you can optimize your fantasy football lineup using mathematical programming.   Full code below:   proc optmodel; set <str> PLAYERS; set IVARS = /'nHits' 'nBB'/; num y{PLAYERS}; num x{PLAYERS,IVARS}; read data sashelp.baseball into PLAYERS=[Name] y=nRuns {k in IVARS} <x[Name,k]=col(k)>; var Intercept; var Beta{IVARS}; min Obj = sum{i in PLAYERS} (y[i] - (Intercept + sum{k in IVARS} Beta[k]*x[i,k]))**2; solve; print Intercept; print Beta; create data work.output from [player] = {i in PLAYERS} y pred=(Intercept + sum{k in IVARS} Beta[k]*x[i,k]); quit;     Find more articles from SAS Global Enablement and Learning here. ... View more

  Okay, let me set the stage for you... You're relaxing on the couch on a lazy fall Sunday afternoon. It's 12:55 pm, and only five short minutes separate you from a full afternoon of NFL football. All is right with the world. As the pregame shows conclude with predictions and final commentary, you come to the horrifying realization that you haven't set your fantasy football lineup yet!   Instantly you spring up from the couch and begin frantically searching for your phone to log in to your fantasy football app, knowing that with each passing second, the window of opportunity to set your starting lineup is closing. Once the clock strikes 1:00 pm, the majority of your roster will be locked. Your fantasy week, and arguably your entire fantasy football season, hangs in the balance on the roster moves you make in the next, now, 3 minutes.   You've managed to locate your phone and fumble your way into the fantasy football app. 2 minutes until 1:00 pm.   On your team's homepage, you see that your current starting lineup is projected to score 88.52 points this week (i.e., your team is represented by the Johnny Cash icon).    Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   Had you have been more proactive, you would've had time to casually check your starting lineup throughout the week to move as many players into or out of the starting lineup as necessary. However, given that you're now in the 11th hour (and 59th minute), you have two options:   Keep your starting roster as-is, with a projected score of 88.52 Select the green "OPTIMIZE" button to automatically set your lineup, which, according to the display, will increase your team's projected score by +8.45 points to 96.97   It's not a trick question, I promise. It's like being asked, Would you rather win an $88,520 jackpot or a $96,970 jackpot? It's a no-brainer.   Without hesitation, you select the big, green "OPTIMIZE" button, save your lineup, breathe a momentous sigh of relief, and hop back on the couch to catch the 1:00 pm kickoff. Problem solved. Crisis averted.     This series will walk through the mechanics of the optimization model embedded within the green "OPTIMIZE" button, allowing you to instantly optimize your fantasy football lineup. Throughout this series, you will learn:   How to formulate and solve this model using different solvers within the OPTMODEL procedure in SAS Optimization (and how to choose the best one!) Common programming syntax used in the OPTMODEL procedure (it's a little bit different than what you might be used to!) Automatic detection and LP relaxation strategies for pure network problems (and, what exactly is a "pure network" problem?!) What reduced costs are in linear programming (and why you should care!) Decomposition techniques for solving large optimization models (after all, Viya is built to solve problems on big data!) How to solve this problem using the OPTNETWORK procedure (and yes, believe it or not, this problem can be formulated as a network optimization model!)   Each subsequent post in this series will introduce or expound upon one or more of these topics in an applied, approachable way. You needn't be an expert to understand and apply the concepts discussed throughout this series. In fact, this series is specifically written for data scientists and quantitative professionals with minimal exposure to mathematical optimization, so consider this your crash course in getting up to speed with the tools and functionality packed into SAS Optimization!   But first thing's first. Let's break down the problem that the big green "OPTIMIZE" button is attempting to solve.   __________________   In my fantasy football league, a roster is composed of 15 players. Each player plays one of the following positions: Quarterback (QB) Running back (RB) Wide receiver (WR) Tight end (TE) Kicker (K) Defense (DEF)   You're allowed to carry as many players at each position as you wish, so long as the total number of players does not exceed 15.   The starting roster is composed of 9 players with the following position limits:   1 QB 2 RBs 2 WRs 1 TE 1 flex position (W/R) containing either a RB or WR 1 K 1 DEF   My current starting lineup, which respects the position limits, is shown below:     Highlighted in yellow to the right of each player is their projected point total for the current week. For example, my starting quarterback, D. Carr, is projected to score 16.53 points this week. One of my running backs, S. Perine, is projected to score 6.43 points this week, and so on.   Summing the projected points across all players in the starting lineup yields a total projected team score of 88.52, which is my team's projected point total for the week previously shown on the homepage above (next to the Johnny Cash icon).   The remaining 6 players are on the bench (i.e., not in the starting lineup). Similar to the players in the starting lineup, each player on the bench has a projected point total for the current week. However, because they're on the bench, their points do not count towards the team's score for the week.     The +8.45 value shown to the left of the green "OPTIMIZE" button indicates that the projected point total for the week will increase by 8.45 points (i.e., from 88.52 currently, to 96.97) simply by clicking that button, which will automatically reallocate the players in my starting lineup to maximize the projected point total for the current week.   The problem is small enough that we can manually optimize the starting lineup. In other words, which players should be swapped to increased the projected point total to 96.97?   Scanning back and forth between the starting lineup and the bench, in no particular order, we can replace:   <RB> S. Perine with <RB> J. Mixon (+7.54) <TE> D. Schultz with <TE> K. Pitts (+0.36) <K> W. Lutz with <K> J. Tucker (+0.55)   Making these three changes will result in a total point gain of +8.45, increasing the team's projected point total to 96.97 for the current week. As you'll learn in a later post, these values (e.g., +7.54, +0.36, +0.55) are defined as the reduced costs in a linear programming model, but let's not get ahead of ourselves.    We just did manually what clicking the green "OPTIMIZE" button will do automatically. However, the reality is there are tens of thousands of fantasy football players spread across numerous leagues. Providing the opportunity to OPTIMIZE each team's lineup every week requires an efficient, automated approach, and mathematical programming is the right tool to perform this task.    To formulate this as a mathematical programming model in SAS Optimization, let's first start with the data.   The data set below lists each position, the maximum allowed at each position in the starting lineup, along with a binary yes/no indicator if that position is eligible for the flex position.   data positions; input pos $ max flex_pos $; datalines; QB 1 no RB 3 yes WR 3 yes TE 1 no K 1 no DEF 1 no ;   The second data set lists each player on the roster, along with their position, projected points, opponent, current health status, and current status as either a starter or on the bench.   data roster; input player $ proj_points pos $ opp $ status $ curr_allocation $; datalines; Mayfield 15.22 QB @NO H bench Mixon 13.97 RB @KC H bench Jacobs 13.43 RB @IND Q starter Sutton 9.43 WR LAC Q starter Waddle 13.13 WR @BAL Q starter Schultz 7.67 TE TEN H starter Kirk 0.00 WR CAR IR bench Carr 16.53 QB @TB H starter Perine 6.43 RB LAC H starter Boyd 7.54 WR @KC H starter Dell 0.00 WR @NYJ IR bench Pitts 8.03 TE @CHI H bench Tucker 7.12 K MIA H bench Lutz 6.57 K LAC H starter SF49 7.79 DEF @WAS H starter ;   The SQL code below modifies the roster data set to flag all RBs and WRs as flex eligible.   proc sql; create table work.roster_wflexelig as select * , case when pos in(select distinct pos from positions where flex_pos='yes') then 'yes' else 'no' end as flex_eligible from roster; drop table roster; quit;   Within the OPTMODEL procedure, the first step declares a set of players, declares the parameters indexed by players, and reads the roster_wflexelig data set into the OPTMODEL procedure.   proc optmodel; set <str> PLAYERS; num proj_points{PLAYERS}; num starter{PLAYERS}; str pos{PLAYERS}; str opp{PLAYERS}; str status{PLAYERS}; str curr_allocation{PLAYERS}; str flex_eligible{PLAYERS}; read data roster_wflexelig into PLAYERS=[player] proj_points pos opp status curr_allocation flex_eligible;   Next, a set of positions is declared, followed by the max and flex_pos parameters indexed by the set of positions. In other words, there is one max value for each position, and one flex_pos for each position, just like the in the data set.   set <str> POSITIONS; num max{POSITIONS}; str flex_pos{POSITIONS}; read data positions into POSITIONS = [pos] max flex_pos;   One modeling approach to solve this problem is to declare a set of binary decision variables. In other words, for each player, is player p assigned to the starting lineup: yes (1) or no (0)?   var Assign{PLAYERS} binary;   The objective function maximizes the total projected points across all players in the starting lineup.   max Total_Points = sum{p in PLAYERS} proj_points[p]*Assign[p];   The first constraint limits the total number players at each position in the starting lineup to be less than or equal to the maximum allowed.   con pos_limit{po in POSITIONS}: sum{p in PLAYERS: pos[p] = po} Assign[p] <= max[po];   To accommodate the RB, WR, and flex positions, there can either be 2 RBs and 3 WRs in the starting lineup, or 3 RBs and 2 WRs. In other words, the total number of RBs and WRs in the starting lineup can't exceed five.   Within the curly brackets below, the colon operator : acts like a where clause, so this constraint is effectively summing across all RBs and WRs, limiting the total number assigned to the starting lineup to be less than or equal to five.   con flex_position: sum{p in PLAYERS: flex_eligible[p] = 'yes'} Assign[p] <= 5;   We'll investigate several decomposition options within the solve statement in subsequent parts of this series, but for this first post, we'll simply call the default solve statement and let the OPTMODEL procedure determine and apply the most appropriate default solver.   solve;   The delta parameter calculates how many more projected points my team is expected to earn compared to my lineup as it's currently set.   num delta = Total_Points.sol - sum{p in PLAYERS: curr_allocation[p] = 'starter'} proj_points[p]; print delta;   Lastly, two output data tables are created containing my optimal starting and bench lineups.   create data optimal_starting_lineup from [player] = {p in PLAYERS: Assign[p] > 0.5} assign curr_allocation pos proj_points opp status; create data optimal_benched_lineup from [player] = {p in PLAYERS: Assign[p] < 0.5} assign curr_allocation pos proj_points opp status; quit;   The Solution Summary table shows the mixed-integer linear programming (MILP) solver was chosen by default, and the projected point total is maximized at 96.97. The delta parameter tells us the optimal starting lineup is projected to yield +8.45 points more than my current starting lineup.     The optimal_starting_lineup table, created from the OPTMODEL procedure, contains the starting lineup yielding the largest possible projected point total of 96.97. Notice the optimal solution makes the same three changes that we made manually above.     The next post in this series will decompose and solve this same problem using Dantzig-Wolfe decomposition. Stay tuned to learn the numerous ways you can use the OPTMODEL procedure to optimize your fantasy football lineup!   Full OPTMODEL code below:   proc optmodel; set <str> PLAYERS; num proj_points{PLAYERS}; num starter{PLAYERS}; str pos{PLAYERS}; str opp{PLAYERS}; str status{PLAYERS}; str curr_allocation{PLAYERS}; str flex_eligible{PLAYERS}; read data roster_wflexelig into PLAYERS=[player] proj_points pos opp status curr_allocation flex_eligible; set <str> POSITIONS; num max{POSITIONS}; str flex_pos{POSITIONS}; read data positions into POSITIONS = [pos] max flex_pos; var Assign{PLAYERS} binary; max Total_Points = sum{p in PLAYERS} proj_points[p]*Assign[p]; con pos_limit{po in POSITIONS}: sum{p in PLAYERS: pos[p] = po} Assign[p] <= max[po]; con flex_position: sum{p in PLAYERS: flex_eligible[p] = 'yes'} Assign[p] <= 5; solve; num delta = Total_Points.sol - sum{p in PLAYERS: curr_allocation[p] = 'starter'} proj_points[p]; print delta; create data optimal_starting_lineup from [player] = {p in PLAYERS: Assign[p] > 0.5} assign curr_allocation pos proj_points opp status; create data optimal_benched_lineup from [player] = {p in PLAYERS: Assign[p] < 0.5} assign curr_allocation pos proj_points opp status; quit;     Find more articles from SAS Global Enablement and Learning here. ... View more

The objective of this article is to introduce the PageRank algorithm, and learn how it can be used to rank local, college, and professional sports teams in a league. This article describes the intuition behind the algorithm, the matrix algebra for computing PageRank scores (i.e., team rankings), and demonstrates two ways to execute the algorithm in SAS Viya.   The PageRank algorithm was originally developed by Google's co-founders Larry Page and Sergey Brin to rank the importance of web pages on the internet. The algorithm is based on the premise that the internet can be represented as one giant network, where the nodes represent web pages, and the directed links represent hyperlinks connecting one web page to another.   For example, on the https://www.mlb.com webpage, at the bottom is a hyperlink connecting to the Boys and Girls Clubs of America (https://www.bgca.org/).     Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   Constructed as a directed network, that's equivalent to an outbound link from https://www.mlb.com (node A) pointing to https://www.bgca.org (node B).    For ranking the importance of nodes in a network, especially for giant networks with numerous inbound and outbound links, one way that nodes could be ranked is to simply count the number of inbound links each node has. The idea is that the more inbound links a node has, the more important it is, and as a result, it should receive a higher rank compared to other nodes in the network.   This ranking method is called degree centrality, and it's a simple, straightforward way to rank node importance. In some networks it can be quite effective. Kim Kardashian currently has over 350 million followers (i.e., over 350 million inbound links) on social media, and it's why companies pay her hundreds of thousands of dollars to promote their product line(s). The more followers she has, the larger of an audience she reaches when posting on social media. More followers = more inbound links = higher in-degree centrality = more important in that network, at least from a marketing standpoint.   PageRank centrality works differently. For the PageRank algorithm, the number of inbound links is only part of the story.   Consider a Winner network, where each node represents a team in a college or professional sports league. Inbound links represent wins, and outbound links represent losses.   Counting the inbound and outbound links for each node simply counts the number of wins and losses for each team, which, if you recall from a moment ago, is equivalent to a team’s in- and out-degree centrality.     Newspapers and sports websites have been displaying standings tables like these forever. However, notice that by constructing the teams, wins, and losses as a network, we get some additional information that the standings tables (i.e., in- and out-degree centrality measures) ignore. Constructed as a network, not only can wins and losses be counted, but the network also tells us who defeated whom. This is critical to understanding how the PageRank algorithm works, because, unlike the standings tables that only tally wins and losses, the PageRank algorithm does not consider all inbound and outbound links to be equally as important to a team's ranking.   For example, defeating a strong opponent (i.e., one with many inbound links, and few outbound links) improves a team's ranking more than defeating a weaker opponent. Similarly, losing to a strong opponent penalizes a team's ranking less than losing to a weaker opponent.   In other words, the PageRank algorithm not only considers the number of in- and outbound links (i.e., wins and losses) each team has, it also factors in the quality of each team's wins and losses when determining its rank.   In the PageRank algorithm, a fundamental tradeoff exists between these two factors (i.e., total wins and losses vs. opponent quality). There's a famous quote by economist, Thomas Sowell, "There are no solutions, there are only trade-offs; and the goal is to find the best trade-off you can." While Sowell was referring to economic policy, it's equally true for the PageRank algorithm.   So, then, how does the PageRank algorithm determine which factor is more important to a team’s ranking? The answer is that it depends on what you (yes, you... this is a user-defined parameter) set the damping parameter, alpha, to be.      The closer alpha is to 1, the more the emphasis is placed on opponent quality. When alpha is set to 0.85, notice below that Team C (1 win, 1 loss, .500 winning percentage) has a higher PageRank value than Team E (6 wins, 3 losses, .667 winning percentage).   The reason is Team C's lone win is against the best team in the league (i.e., the one with the most inbound links (wins) and fewest outbound links (losses)), Team B. Even though Team E has a higher winning percentage, the overwhelming majority of Team E's wins are against inferior opponents.   As the damping parameter, alpha, approaches 0, more emphasis is placed on total inbound links (i.e., wins), and less on opponent quality. For example, when alpha is changed to 0.25, Team E is now ranked higher than Team C, simply because Team E has more overall wins (i.e., inbound links), regardless of opponent quality.   In its original context, alpha, represented the probability of traversing the internet network via hyperlink, and (1 - alpha) represented the probability of "teleporting" from one node to the next. Think about it: we don't always traverse the internet by hopping from hyperlink to hyperlink. There's nothing preventing us from going to the search bar, typing in our desired web page address, and "teleporting" across the network to the website of our choice.   The point here is twofold. First, the PageRank algorithm can be used on a variety of different types of networks to rank node importance. Second, depending on what the nodes and links in the network represent, the damping parameter can take on different interpretations, but the underlying mechanics for ranking nodes using the PageRank algorithm remain consistent.   Ranking nodes using the PageRank algorithm has never been easier in SAS Viya. The NETWORK procedure includes a number of graph theory and network analysis algorithms to augment data mining and machine learning approaches. PageRank is one of the many network algorithms embedded within the NETWORK procedure, and it's easy to use.   The code below reads in a directed edge list, computes PageRank centrality, and outputs the results to a new CAS table called PageRankScores. The pageRankalpha option allows you to set the damping parameter to a user-defined level between 0 and 1.   proc network direction = directed links = mycas.linksTable outNodes = mycas.PageRankScores; centrality pageRank = unweight pageRankalpha = 0.85; run;   Below are the results from the CAS table defined in the outNodes= statement. Nodes are sorted by most important node. While not shown below, sometimes it’s preferred to normalize the PageRank scores by dividing each score by the largest value, so the node with the largest PageRank score equals 1, and all other nodes receive a proportion of its value.     Similar to other AI and machine learning algorithms, behind the scenes, the PageRank algorithm is computed using linear algebra. SAS has a rich and comprehensive matrix programming language called SAS IML (Interactive Matrix Language), allowing users to build and customize AI and ML applications.   Within the iml Action in SAS Viya, the code below uses the methodology proposed in Langville and Meyer, 2012 to construct a custom pageRank( ) function to compute PageRank scores from an adjacency matrix.   start pageRank(AdjMatrix, alpha=0.85, trProbVec=); n = nrow(AdjMatrix); if isskipped(trProbVec) then do; e = j(n,1,1); eT = t(e); vT = (1/n)*eT; customTeleVec = 'No'; end; else if ^isskipped(trProbVec) then do; vT = trProbVec; customTeleVec = 'Yes'; end; H = AdjMatrix / (j(n,1,1e-12) <> AdjMatrix[,+]); leafBinary=j(n,1,0); /* set leafBinary[i] = 0 for all i */ if ^isempty(loc(AdjMatrix[,+]=0)) then leafBinary[loc(AdjMatrix[,+]=0)]=1; pi = vT; /* initialize page rank vector to equal vT */ maxIters = 1e8; tolerance = 1e-12; i = 0; /* initialize iteration scalar */ t = time(); do while (i < maxIters); /* power method */ i = i + 1; prevpi = pi; pi = alpha * pi * H + (alpha * (pi * leafBinary) + 1 - alpha) * vT; delta = distance(pi, prevpi) / max(prevpi[+,]); if delta < tolerance then do; iters = i; i = maxIters; pageRankVec = t(pi); end; end; time = time() - t; /* calculate total time */ return(pageRankVec); finish;   Given an adjacency matrix, A, representing the nodes and links in a network, a new matrix, H, is created by converting the elements (i.e., links) into transition probabilities by dividing by the row sum (i.e., total outbound links) for each node. The result is a new, row stochastic matrix, H, where each row sums to 1. But wait a second... the first row in H doesn't sum to 1! It turns out that in this network, the node corresponding to row 1 in the adjacency matrix doesn't have any outbound links, so dividing each element in row 1 by the row sum is equivalent to dividing 0 by 0. This node is called a dangling node. In the Winner network above, node A is an example of a dangling node because it has no outbound links (i.e., losses).   In a connected network, dangling nodes teleport with equal probability (1/n*eT) to each node in the network, where eT is a row vector of ones and n is the total number of nodes in the network. In other words, for dangling nodes, every element receives the value (1/n*eT), meaning that each dangling node is equally likely to teleport to any other node in the network.   To accommodate dangling nodes, a new matrix, S, is created, where S = H + a(1/n * eT). The elements in the nx1 binary vector, a, equal one if the node is a dangling node, and zero otherwise. Therefore, a(1/n * eT) is applied only to the dangling nodes in the network. If the network has no dangling nodes, the equation effectively reduces to S = H.     Working through the equation, now every row in S sums to 1 regardless of its dangling node status.   The damping parameter, alpha, is introduced by converting the S matrix to the G matrix, where G = alpha * S + (1 - alpha) * 1/n * e * eT. The alpha value used to create the G matrix below is 0.85.     Notice by converting S to G, transition probabilities are dampened based on the value of alpha, so that all transition probabilities now contain non-zero values.     The PageRank scores are the dominant left-hand eigenvector of the G matrix, but because G is completely dense, and H is sparse, it’s computationally advantageous to solve G in terms of H, where G = alpha * H + (alpha * a + (1 - alpha) * e) * 1/n * eT. The SAS IML function created above uses the power method, which is an iterative technique, to find the dominant eigenpair from the matrix, and is one of several options available within the NETWORK procedure.   Inside of the IML action, the adjacency matrix from the Winner network is created, along with a vector of row and column names, and the pageRank() function is used to compute PageRank scores.    nodeVec = {'A' 'B' 'C' 'D' 'E' 'F' 'G' 'H' 'I' 'J' 'K'}; A = { 0 0 0 0 0 0 0 0 0 0 0, 0 0 1 0 0 0 0 0 0 0 0, 0 1 0 0 0 0 0 0 0 0 0, 1 1 0 0 0 0 0 0 0 0 0, 0 1 0 1 0 1 0 0 0 0 0, 0 1 0 0 1 0 0 0 0 0 0, 0 0 0 0 1 0 0 0 0 0 0, 0 0 0 0 1 0 0 0 0 0 0, 0 1 0 0 1 0 0 0 0 0 0, 0 1 0 0 1 0 0 0 0 0 0, 0 1 0 0 1 0 0 0 0 0 0 }; PageRankVec = pageRank(A,0.85); PageRankVecNorm = PageRankVec / PageRankVec[<>]; print PageRankVec[rowname=nodeVec] PageRankVecNorm[rowname=nodeVec];   The output from the print statement below shows consistent results between the two approaches in SAS Viya.   Interested in learning more about the functionality packed into the NETWORK procedure, like community detection, shortest paths, network projection, pattern matching, cycles, etc., all with real world applications? Check out our Network Analysis and Network Optimization in SAS Viya training course.   See the full PROC NETWORK and SAS IML programs below.   SAS IML:     cas mysess; libname mycas cas sessref=mysess; proc cas; session mysess; loadactionset 'iml'; /* load the IML action set */ source PageRank; start pageRank(AdjMatrix, alpha=0.85, trProbVec=); n = nrow(AdjMatrix); if isskipped(trProbVec) then do; e = j(n,1,1); eT = t(e); vT = (1/n)*eT; customTeleVec = 'No'; end; else if ^isskipped(trProbVec) then do; vT = trProbVec; customTeleVec = 'Yes'; end; H = AdjMatrix / (j(n,1,1e-12) <> AdjMatrix[,+]); leafBinary=j(n,1,0); /* set leafBinary[i] = 0 for all i */ if ^isempty(loc(AdjMatrix[,+]=0)) then leafBinary[loc(AdjMatrix[,+]=0)]=1; pi = vT; /* initialize page rank vector to equal vT */ maxIters = 1e8; tolerance = 1e-12; i = 0; /* initialize iteration scalar */ t = time(); do while (i < maxIters); /* power method */ i = i + 1; prevpi = pi; pi = alpha * pi * H + (alpha * (pi * leafBinary) + 1 - alpha) * vT; delta = distance(pi, prevpi) / max(prevpi[+,]); if delta < tolerance then do; iters = i; i = maxIters; pageRankVec = t(pi); end; end; time = time() - t; /* calculate total time */ print iters time; return(pageRankVec); finish; nodeVec={'A' 'B' 'C' 'D' 'E' 'F' 'G' 'H' 'I' 'J' 'K'}; A = { 0 0 0 0 0 0 0 0 0 0 0, 0 0 1 0 0 0 0 0 0 0 0, 0 1 0 0 0 0 0 0 0 0 0, 1 1 0 0 0 0 0 0 0 0 0, 0 1 0 1 0 1 0 0 0 0 0, 0 1 0 0 1 0 0 0 0 0 0, 0 0 0 0 1 0 0 0 0 0 0, 0 0 0 0 1 0 0 0 0 0 0, 0 1 0 0 1 0 0 0 0 0 0, 0 1 0 0 1 0 0 0 0 0 0, 0 1 0 0 1 0 0 0 0 0 0 }; PageRankVec = pageRank(A,0.85); PageRankVecNorm = PageRankVec / PageRankVec[<>]; print PageRankVec[rowname=nodeVec] PageRankVecNorm[rowname=nodeVec]; endsource; iml / code=PageRank nthreads=4; run;   PROC NETWORK:   cas mysess; libname mycas cas sessref=mysess; data mycas.linksTable; input from $ to $; datalines; B C C B D A D B E B E D E F F B F E G E H E I E I B J E J B K B K E ; proc network direction = directed links = mycas.linksTable outNodes = mycas.PageRankScores; centrality pageRank = unweight pageRankalpha = 0.85; run;   Find more articles from SAS Global Enablement and Learning here. ... View more