elvis/cmdla
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

first commit with current lessons

Browse Source
This commit is contained in:
elvis
2023-10-29 02:06:02 +01:00
parent ececf6cfcb
commit b3fc9cb021
204 changed files with 31603 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1224
10-04/GMQ.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

247
10-04/GMQ.jl Normal file
View File

@ -0,0 +1,247 @@
using LinearAlgebra
using Printf
using Plots
function gmq(Q::Matrix, q::Vector; x::Union{Vector, Nothing}=nothing, fStar::Real=-Inf, alpha::Real=0 , MaxIter::Int=1000 , eps::Real=1e-6, plt::Union{Plots.Plot, Nothing}=nothing, Plotf::Int=2, printing::Bool=true)::Tuple{Vector, String}
# Plotf
# 0 = nothing is plotted
# 1 = the function value / gap are plotted
# 2 = the level sets of f and the trajectory are plotted (when n = 2)
Interactive = true # if we pause at every iteration
Streamlined = true # if the streamlined version of the algorithm, with
# only one O( n^2 ) operation per iteration, is used
# reading and checking input- - - - - - - - - - - - - - - - - - - - - - - -
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if !isreal(Q)
throw(ArgumentError(Q, "Q not a real matrix"))
end
n = size(Q, 1)
if n <= 1
throw(ArgumentError(Q, "Q is too small"))
end
if n != size(Q, 2)
throw(ArgumentError(Q, "Q is not square"))
end
if !isreal(q)
throw(ArgumentError(q, "q not a real vector"))
end
if size(q, 1) != n
throw(ArgumentError(q, "q size does not match with Q"))
end
if x == nothing
x = zeros(n, 1)
end
if !isreal(x)
throw(ArgumentError(x, "x not a real vector"))
end
if size(x, 1) != n
throw(ArgumentError(x, "x size does not match with Q"))
end
if MaxIter < 1
throw(ArgumentError(MaxIter, "MaxIter too small"))
end
if eps < 0
throw(ArgumentError(eps, "eps can not be negative"))
end
# initializations - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if printing
print("Gradient method for quadratic functions ")
if alpha == 0
print("(optimal stepsize)\n")
else
print("(fixed stepsize)\n")
end
print("iter\tf(x)\t\t\t||g||")
end
if fStar > - Inf
if printing
print("\t\tgap\t\trate")
end
prevf = Inf
end
if printing
if alpha == 0
print("\t\talpha")
end
print("\n\n")
end
i = 0;
if Plotf == 1
gap = []
end
if Streamlined
g = Q * x + q
end
if Plotf == 1 && plt == nothing
plt = plot(yscale = :log,
xlims=(0, MaxIter),
ylims=(1e-15, Inf),
guidefontsize=16)
elseif Plotf == 2 && plt == nothing
plt = plot()
end
status = ""
# main loop - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
while true
if !Streamlined
g = Q * x + q
end
ng = norm(g)
f = dot((g + q)', x) / 2 # 1/2 x^T Q x + q x
# = 1/2 ( x^T Q x + 2 q x )
# = 1/2 x^T ( Q x + q + q )
# = 1/2 ( q + g ) x
i += 1
if printing
@printf("%4d\t%1.8e\t\t%1.4e", i, f, ng)
end
if fStar > -Inf
gapk = (f - fStar)/maximum([abs(fStar), 1])
if printing
@printf("\t%1.4e", gapk)
if prevf < Inf
@printf("\t%1.4e", (f - fStar)/(prevf - fStar))
else
@printf("\t\t")
end
end
prevf = f
if Plotf == 1
push!(gap, gapk)
end
end
# stopping criteria - - - - - - - - - - - - - - - - - - - - - - - - - -
if ng <= eps
status = "optimal"
if alpha == 0 && printing
print("\n")
end
break
end
if i > MaxIter
status = "stopped"
if alpha == 0 && printing
print("\n")
end
break
end
# compute step size - - - - - - - - - - - - - - - - - - - - - - - - - -
# meanwhile, check if f is unbounded below
# note that if alpha > 0 this is only used for the unboundedness check
# which is a bit of a waste, but there you go; anyway, in the
# streamlined version this only costs O( n )
if Streamlined
v = Q * g;
den = dot(g', v)
else
den = dot(g', Q * g)
end
if den <= 1e-14
# this is actually two different cases:
# - g' * Q * g = 0, i.e., f is linear along g, and since the
# gradient is not zero, it is unbounded below
#
# - g' * Q * g < 0, i.e., g is a direction of negative curvature for
# f, which is then necessarily unbounded below
if printing
if alpha == 0
print("\n")
end
@printf("g' * Q * g = %1.4e ==> unbounded\n", den)
end
status = "unbounded"
break
end
if alpha > 0
t = alpha
else
t = ng^2 / den # stepsize
if printing
@printf("\t%1.2e", t)
end
end
if printing
print("\n")
end
# compute new point - - - - - - - - - - - - - - - - - - - - - - - - - -
# possibly plot the trajectory
if n == 2 && Plotf == 2
PXY = hcat(vec(x), vec(x - t * g))
plot!(PXY[1,:],
PXY[2,:],
linestyle=:solid,
linewidth=2,
markershape=:circle,
seriescolor=colorant"black",
label="")
end
x = x - t * g
if Streamlined
g = g - t * v
end
if Interactive
#readline()
end
if Plotf != 0
#IJulia.clear_output(true)
#display(plt)
end
end
if Plotf == 1
plot!(plt,
gap,
linewidth=2,
seriescolor=colorant"black")
display(plt)
elseif Plotf == 2
display(plt)
end
(vec(x), status)
end

319
10-04/GMQ.m Normal file
View File

@ -0,0 +1,319 @@
function [ x , status ] = GMQ( Q , q , varargin )
%function [ x , status ] = GMQ( Q , q , x , fStar , alpha , MaxIter , eps )
%
% Apply the Gradient Method (a.k.a., Steepest Descent algorithm) to the
% minimization of the quadratic function
%
% f( x ) = 1/2 x^T Q x + q x
%
% Input:
%
% - Q ([ n x n ] real symmetric matrix, not necessarily positive
% semidefinite): the quadratic part of f
%
% - q ([ n x 1 ] real column vector): the linear part of f
%
% - x ([ n x 1 ] real column vector or empty, optional): the point where to
% start the algorithm from; if not provided or empty, the all-0 n-vector
% is used
%
% - fStar (real scalar, optional, default value Inf): optimal value of f.
% if a non-Inf value is provided it is used to print out stasistics about
% the convergence speed of the algorithm
%
% - alpha (real scalar, optional, default value 0): if alpha > 0, then the
% fixed-stepsize version of the algorithm is run with alpha as the fixed
% stepsize, otherwise the standard exact line search is used
%
% - MaxIter (integer scalar, optional, default value 1000): the maximum
% number of iterations
%
% - eps (real scalar, optional, default value 1e-6): the accuracy in the
% stopping criterion: the algorithm is stopped when the norm of the
% gradient is less than or equal to eps
%
% Output:
%
% - x ([ n x 1 ] real column vector): either the best solution found so far
% (possibly the optimal one) or a direction proving the problem is
% unbounded below, depending on case
%
% - status (string): a string describing the status of the algorithm at
% termination
%
% = 'optimal': the algorithm terminated having proven that x is a(n
% approximately) optimal solution, i.e., the norm of the gradient at x
% is less than the required threshold
%
% = 'unbounded': the algorithm terminated having proven that the problem
% is unbounded below: x contains a direction along which f is
% decreasing to - Inf, either because f is linear along x and the
% directional derivative is not zero, or because x is a direction with
% negative curvature
%
% = 'stopped': the algorithm terminated having exhausted the maximum
% number of iterations: x is the best solution found so far, but not
% necessarily the optimal one
%
%{
=======================================
Author: Antonio Frangioni
Date: 26-12-22
Version 0.2
Copyright Antonio Frangioni
=======================================
%}
Plotf = 2;
% 0 = nothing is plotted
% 1 = the function value / gap are plotted
% 2 = the level sets of f and the trajectory are plotted (when n = 2)
Interactive = true; % if we pause at every iteration
Streamlined = true; % if the streamlined version of the algorithm, with
% only one O( n^2 ) operation per iteration, is used
% reading and checking input- - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if ~ isreal( Q )
error( 'Q not a real matrix' );
end
n = size( Q , 1 );
if n <= 1
error( 'Q is too small' );
end
if n ~= size( Q , 2 )
error( 'Q is not square' );
end
if ~ isreal( q )
error( 'q not a real vector' );
end
if size( q , 2 ) ~= 1
error( 'q is not a (column) vector' );
end
if size( q , 1 ) ~= n
error( 'q size does not match with Q' );
end
if isempty( varargin ) || isempty( varargin{ 1 } )
x = zeros( n , 1 );
else
x = varargin{ 1 };
if ~ isreal( x )
error( 'x not a real vector' );
end
if size( x , 2 ) ~= 1
error( 'x is not a (column) vector' );
end
if size( x , 1 ) ~= n
error( 'x size does not match with Q' );
end
end
if length( varargin ) > 1
fStar = varargin{ 2 };
if ~ isreal( fStar ) || ~ isscalar( fStar )
error( 'fStar is not a real scalar' );
end
else
fStar = - Inf;
end
if length( varargin ) > 2
alpha = varargin{ 3 };
if ~ isreal( alpha ) || ~ isscalar( alpha )
error( 'alpha is not a real scalar' );
end
else
alpha = 0;
end
if length( varargin ) > 3
MaxIter = round( varargin{ 4 } );
if ~ isscalar( MaxIter )
error( 'MaxIter is not an integer scalar' );
end
if MaxIter < 1
error( 'MaxIter too small' );
end
else
MaxIter = 1000;
end
if length( varargin ) > 4
eps = varargin{ 5 };
if ~ isreal( eps ) || ~ isscalar( eps )
error( 'eps is not a real scalar' );
end
if eps < 0
error( 'eps can not be negative' );
end
else
eps = 1e-6;
end
% initializations - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
fprintf( 'Gradient method for quadratic functions ' );
if alpha == 0
fprintf( '(optimal stepsize)\n' );
else
fprintf( '(fixed stepsize)\n' );
end
fprintf( 'iter\tf(x)\t\t\t||g||');
if fStar > - Inf
fprintf( '\t\tgap\t\trate');
prevf = Inf;
end
if alpha == 0
fprintf( '\t\talpha' );
end
fprintf( '\n\n' );
i = 0;
if Plotf == 1
gap = [];
end
if Streamlined
g = Q * x + q;
end
% main loop - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
while true
% compute function value and gradient - - - - - - - - - - - - - - - - -
if ~ Streamlined
g = Q * x + q;
end
ng = norm( g );
f = ( g + q )' * x / 2; % 1/2 x^T Q x + q x = 1/2 ( x^T Q x + 2 q x )
% = 1/2 x^T ( Q x + q + q ) = 1/2 ( q + g ) x
i = i + 1;
% output statistics - - - - - - - - - - - - - - - - - - - - - - - - - - -
fprintf( '%4d\t%1.8e\t\t%1.4e' , i , f , ng );
if fStar > - Inf
gapk = ( f - fStar ) / max( [ abs( fStar ) 1 ] );
fprintf( '\t%1.4e' , gapk );
if prevf < Inf
fprintf( '\t%1.4e' , ( f - fStar ) / ( prevf - fStar ) );
else
fprintf( '\t\t' );
end
prevf = f;
if Plotf == 1
gap( end + 1 ) = gapk;
semilogy( gap , 'Color' , 'k' , 'LineWidth' , 2 );
xlim( [ 0 MaxIter ] );
ylim( [ 1e-15 inf ] );
ax = gca;
ax.FontSize = 16;
ax.Position = [ 0.03 0.07 0.95 0.92 ];
ax.Toolbar.Visible = 'off';
end
end
% stopping criteria - - - - - - - - - - - - - - - - - - - - - - - - - -
if ng <= eps
status = 'optimal';
if alpha == 0
fprintf( '\n' );
end
break;
end
if i > MaxIter
status = 'stopped';
if alpha == 0
fprintf( '\n' );
end
break;
end
% compute step size - - - - - - - - - - - - - - - - - - - - - - - - - -
% meanwhile, check if f is unbounded below
% note that if alpha > 0 this is only used for the unboundedness check
% which is a bit of a waste, but there you go; anyway, in the
% streamlined version this only costs O( n )
if Streamlined
v = Q * g;
den = g' * v;
else
den = g' * Q * g;
end
if den <= 1e-14
% this is actually two different cases:
%
% - g' * Q * g = 0, i.e., f is linear along g, and since the
% gradient is not zero, it is unbounded below
%
% - g' * Q * g < 0, i.e., g is a direction of negative curvature for
% f, which is then necessarily unbounded below
%
if alpha == 0
fprintf( '\n' );
end
fprintf( 'g'' * Q * g = %1.4e ==> unbounded\n' , den );
status = 'unbounded';
break;
end
if alpha > 0
t = alpha;
else
t = ng^2 / den; % stepsize
fprintf( '\t%1.2e' , t );
end
fprintf( '\n' );
% compute new point - - - - - - - - - - - - - - - - - - - - - - - - - -
% possibly plot the trajectory
if n == 2 && Plotf == 2
PXY = [ x , x - t * g ];
line( 'XData' , PXY( 1 , : ) , 'YData' , PXY( 2 , : ) , ...
'LineStyle' , '-' , 'LineWidth' , 2 , 'Marker' , 'o' , ...
'Color' , [ 0 0 0 ] );
end
x = x - t * g;
if Streamlined
g = g - t * v;
end
% iterate - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if Interactive
pause;
end
end
% end of main loop- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end % the end- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

111
10-04/genQf.jl Normal file
View File

@ -0,0 +1,111 @@
using Random
using LinearAlgebra
function genQf(n::Integer; fseed::Integer=0, rank::Real=1.1, conv::Real=1, ecc::Real=0.99, dom::Real=1, box::Real=1, q::Union{Vector, Nothing}=nothing, v::Union{Real, Nothing}=nothing)::Tuple{Matrix, Union{Vector, Nothing}, Union{Real, Nothing}}
if n <= 0
throw(ArgumentError(n, "n must be > 0"))
end
if rank <= 0
throw(ArgumentError(rank, "rank must be > 0"))
end
if !(0 <= conv <= 1)
throw(ArgumentError(conv, "conv must be in [0, 1]"))
end
if !(0 <= ecc < 1)
throw(ArgumentError(ecc, "ecc must be in [0, 1)"))
end
if dom < 0
throw(ArgumentError(dom, "dom must be >= 0"))
end
if box == 0
throw(ArgumentError(box, "box must not be 0"))
end
Random.seed!(fseed)
# generate Q- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# first step: generate the appropriate rank and positive / negative
# definiteness
r = round(Int, rank * n)
p = round(Int, r * conv)
if p > 0
G = rand(p, n)
Q = G' * G
else
Q = zeros(n, n)
end
if r > p
G = rand(r - p, n)
Q = Q - (G' * G)
end
# second step: if dom ~= 1, modify the diagonal
# increase or decrease randomly each element by [ - 1/3 , 1/3 ] of its
# initial value, then multiply it by dom
if dom != 1
D = diag(Q)
D = D .* (1 .+ (2 .* rand(n, 1) .- 1) ./ 3)
D = dom * D
@view(Q[diagind(Q)]) .= D
end
# compute eigenvalue decomposition
F = eigen(Q); # V * D * V' = Q , D( 1 , 1 ) = \lambda_n
V = F.vectors
D = F.values
if D[1] > 1e-14
# modify eccentricity only if \lambda_n > 0, for when \lambda_n = 0 the
# eccentricity is 1 by default
#
# the formula is:
#
# \lambda_i - \lambda_n 2 * ecc
# \lambda_i = \lambda_n + --------------------- * \lambda_n -------
# \lambda_1 - \lambda_n 1 - ecc
#
# This leaves \lambda_n unchanged, and modifies all the other ones
# proportionally so that
#
# \lambda_1 - \lambda_n
# --------------------- = ecc (exercise: check)
# \lambda_1 - \lambda_n
l = D[1] .* ones(n) + ((D[1] / (D[n] - D[1])) * (2 * ecc / (1 - ecc))) .* (D - D[1] .* ones(n))
Q = V * diagm(l) * V';
end
if q != nothing || v != nothing
# if so required generate q- - - - - - - - - - - - - - - - - - - - - - -
#
# we first generate the unconstrained minimum x_* of the problem in the
# box [ - abs( box ) , abs( box ) ] and then we set q = - Q * x_*
x = 2 * abs(box) .* rand(n) .- abs(box)
q = -Q * x
# if so required, we now randomly destroy the alignment between q and
# the image of Q so as to make it hard to solve Q x = q
if ( box < 0 ) && (D[1] <= 1e-14 )
q = q .* ((4/3) .* rand(n) .- (2/3))
end
end
if v != nothing
# if so required compute v. - - - - - - - - - - - - - - - - - - - - -
# v is finite-valued only if either Q is strictly positive definite
# or it is positive semidefinite but q has been constructed in such
# a was that Q * x + q = 0 has a solution (that is the x we have)
if (D[1] > 1e-14) || ((D[1] > -1e-14) && (box > 0))
v = dot(q', x) / 2
else
v = -Inf
end
end
return (Q, q, v)
end

278
10-04/getQf.m Normal file
View File

@ -0,0 +1,278 @@
function [ Q , varargout ] = genQf( n , varargin )
%function [ Q , q , v ] = genQf( n , seed , rank , conv , ecc , dom , box )
%
% Produces the data structure encoding a Quadratic function
%
% f( x ) = (1/2) x^T * Q * x + q * x
%
% Input:
%
% - n (integer, scalar): the size of the problem
%
% - seed (integer, default 0): the seed for the random number generator
%
% - rank (real, scalar, default 1.1): Q will be obtained as Q = G^T G, with
% G a m \times n random matrix with m = rank * n (almost, see "conv"
% below). If rank > 1 then Q can be expected to be full-rank, if
% rank < 1 it will not
%
% - conv (real, scalar, default 1): a parameter controlling whether Q is
% positive (semi)definite, negative (semi)definite, or indefinite. This
% is done by actually computing Q as Q = G_+^T G_+ - G_-^T G_-. The m
% rows of G are partitioned according to conv, i.e., m * conv rows are
% put in G_+ while m * ( 1 - conv ) rows are put in G_-. Hence, conv = 1
% means that Q is positive (semi)definite, conv = 1 means that Q is
% negative (semi)definite, and any value in the middle means that Q is
% indefinite, with "close to 1" meaning "almost positive (semi)definite",
% and "close to 0" meaning "almost negative (semi)definite"
%
% - ecc (real, scalar, default 0.99): the eccentricity of Q, i.e., the
% ratio ( \lambda_1 - \lambda_n ) / ( \lambda_1 + \lambda_n ), with
% \lambda_1 the largest eigenvalue and \lambda_n the smallest one. Note
% that this makes sense only if \lambda_n > 0, for otherwise the
% eccentricity is always 1; hence, this setting is ignored if
% \lambda_n = 0, i.e., Q is not full-rank (see above). An eccentricity of
% 0 means that all eigenvalues are equal, as eccentricity -> 1 the
% largest eigenvalue gets larger and larger w.r.t. the smallest one
%
% - dom (real, scalar, default 1): possibly modifies the matrix by making
% it more or less diagonally dominant. This is obtained by multiplying
% each diagonal element by a random number uniformly picked in the
% interval [ 2 * diag / 3 , 4 * diag / 3 ] (and therefore with average
% value diag). If diag >> 1 this increases the diagonal dominance, if
% diag << 1 this decreases it. The operation is only performed if
% diag ~= 1.
%
% - box (real, scalar, default 1): this parameter controls the generation
% of q. This is done by first randomly generating one (tentatively)
% optimal solution x_* to the unconstrained minimization of f( x ) in the
% symmetric hypercube of size 2 * abs( box ), i.e., by uniformly drawing
% each of its entries in the interval [ - abs( box ) , abs( box ) ];
% then, q is generated as - Q * x_*, so that Q * x_* + q = 0. This
% ensures that f( x ) is bounded below if Q is positive semidefinite, and
% bounded above if Q is negative semidefinite. However, we also want to
% be able to create problems that are positive semidefinite but not
% bounded below, which means that the system Q x + q = 0 must not have a
% solution. This first of all requires Q to be low-rank (some of the
% eigenvalues actually = 0), which can be obtained by putting rank < 1.
% Then, if box < 0 the vector q obtained as above is modified by
% multiplying each of its entries by a random number uniformly picked in
% the interval [ 2 / 3 , 4 / 3 ], which makes it unlikely that the system
% Q x = q still has a solution. However, this makes sense only if Q is
% rank-deficient, hence this is done only if Q is not strictly positive
% definite (for otherwise the system always has a solution).
%
% Output:
%
% - Q: n \times n symmetric real matrix, which is either positive
% (semi)definite, negative (semi)definite, or indefinite according
% to the value of the conv parameter
%
% - q: n \times 1 real vector, optional
%
% - v: real, optional: the optimal value of the optimizattion problem
% min{ f( x ) : x \in \R^n }. This is -\infty if Q is not positive
% semidefinite (see conv above). It is -\infty even if Q is positive
% semidefinite but rank deficient (see rank above) and the system
% Q x + q = 0 has no solution. If instead Q is positive semidefinite
% (whatever its rank) and the system Q x + q = 0 has a solution x_*,
% then x_* is an optimal solution to the problem and
%
% v = f( x_* ) = (1/2) x_*^T * Q * x_* + q * x_*
% = (1/2) [ x_*^T * Q * x_* + q * x_* ] + (1/2) q * x_*
% = (1/2) x_*^T ( Q * x_* + q * ) + (1/2) q * x_*
% = (1/2) q * x_*
%
%{
=======================================
Author: Antonio Frangioni
Date: 15-08-22
Version 0.10
Copyright Antonio Frangioni
=======================================
%}
if ~ isscalar( n ) || ~ isreal( n )
error( 'n not a real scalar' );
end
n = round( n );
if n <= 0
error( 'n must be > 0' );
end
if isempty( varargin )
seed = 0;
else
seed = varargin{ 1 };
if ~ isscalar( seed ) || ~ isreal( seed )
error( 'actv not a real scalar' );
end
seed = round( seed );
end
if length( varargin ) > 1
rank = varargin{ 2 };
if ~ isscalar( rank ) || ~ isreal( rank )
error( 'rank not a real scalar' );
end
if rank <= 0
error( 'rank must be > 0' );
end
else
rank = 1.1;
end
if length( varargin ) > 2
conv = varargin{ 3 };
if ~ isscalar( conv ) || ~ isreal( conv )
error( 'conv not a real scalar' );
end
if ( conv < 0 ) || ( conv > 1 )
error( 'conv must be in [0, 1)' );
end
else
conv = 1;
end
if length( varargin ) > 3
ecc = varargin{ 4 };
if ~ isscalar( ecc ) || ~ isreal( ecc )
error( 'ecc not a real scalar' );
end
if ecc < 0 || ecc >= 1
error( 'ecc must be in [0, 1)' );
end
else
ecc = 0.99;
end
if length( varargin ) > 4
dom = varargin{ 5 };
if ~ isscalar( dom )
error( 'dom not a scalar' );
end
if dom < 0
error( 'dom must be >= 0' );
end
else
dom = 1;
end
if length( varargin ) > 5
box = varargin{ 6 };
if ~ isscalar( box )
error( 'box not a scalar' );
end
if box == 0
error( 'box must not be 0' );
end
else
box = 1;
end
rng( seed );
% generate Q- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% first step: generate the appropriate rank and positive / negative
% definiteness
r = round( rank * n );
p = round( r * conv );
if p > 0
G = rand( p , n );
Q = G' * G;
else
Q = zeros( n , n );
end
if r > p
G = rand( r - p , n );
Q = Q - G' * G;
end
% second step: if dom ~= 1, modify the diagonal
% increase or decrease randomly each element by [ - 1/3 , 1/3 ] of its
% initial value, then multiply it by dom
if dom ~= 1
D = diag( Q );
D = D .* ( 1 + ( 2 * rand( n , 1 ) - 1 ) / 3 );
D = dom * D;
Q = spdiags( D , 0 , Q );
Q = full( Q );
end
% compute eigenvalue decomposition
[ V , D ] = eig( Q ); % V * D * V' = Q , D( 1 , 1 ) = \lambda_n
if D( 1 , 1 ) > 1e-14
% modify eccentricity only if \lambda_n > 0, for when \lambda_n = 0 the
% eccentricity is 1 by default
%
% the formula is:
%
% \lambda_i - \lambda_n 2 * ecc
% \lambda_i = \lambda_n + --------------------- * \lambda_n -------
% \lambda_1 - \lambda_n 1 - ecc
%
% This leaves \lambda_n unchanged, and modifies all the other ones
% proportionally so that
%
% \lambda_1 - \lambda_n
% --------------------- = ecc (exercise: check)
% \lambda_1 - \lambda_n
d = diag( D );
l = d( 1 ) * ones( n , 1 ) + ( d( 1 ) / ( d( n ) - d( 1 ) ) ) * ...
( 2 * ecc / ( 1 - ecc ) ) * ...
( d - d( 1 ) * ones( n , 1 ) );
Q = V * diag( l ) * V';
end
if nargout > 0
% if so required generate q- - - - - - - - - - - - - - - - - - - - - - -
%
% we first generate the unconstrained minimum x_* of the problem in the
% box [ - abs( box ) , abs( box ) ] and then we set q = - Q * x_*
x = 2 * abs( box ) * rand( n , 1 ) - abs( box );
q = - Q * x;
% if so required, we now randomly destroy the alignment between q and
% the image of Q so as to make it hard to solve Q x = q
if ( box < 0 ) && ( D( 1 , 1 ) <= 1e-14 )
q = q .* ( ( 4 / 3 ) * rand( n , 1 ) - 2 / 3 );
end
varargout{ 1 } = q;
if nargout > 1
% if so required compute v. - - - - - - - - - - - - - - - - - - - - -
% v is finite-valued only if either Q is strictly positive definite
% or it is positive semidefinite but q has been constructed in such
% a was that Q * x + q = 0 has a solution (that is the x we have)
if ( D( 1 , 1 ) > 1e-14 ) || ...
( ( D( 1 , 1 ) > -1e-14 ) && ( box > 0 ) )
v = q' * x / 2;
else
v = - Inf;
end
varargout{ 2 } = v;
end
end
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end

383
10-04/lesson.ipynb Normal file
View File

@ -0,0 +1,383 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "a88398a0-03e1-458f-ba48-518fb151b1b3",
"metadata": {},
"outputs": [],
"source": [
"using Plots\n",
"using LinearAlgebra"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "deac13e6-b3dd-4c8d-9c15-504052e60c6d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Matrix{Float64}:\n",
" 0.693988 0.0373642 0.306196\n",
" 0.475249 0.0523086 0.525352\n",
" 0.275428 0.0212671 0.431897"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = rand(3,3)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "a7cfa41c-b406-475b-8a8e-bc56cab68cfa",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}\n",
"U factor:\n",
"3×3 Matrix{Float64}:\n",
" -0.649854 0.735055 0.193349\n",
" -0.624454 -0.371326 -0.687149\n",
" -0.433297 -0.567284 0.700316\n",
"singular values:\n",
"3-element Vector{Float64}:\n",
" 1.1253009779161582\n",
" 0.27878010490969024\n",
" 0.013851237498536341\n",
"Vt factor:\n",
"3×3 Matrix{Float64}:\n",
" -0.770554 -0.0587937 -0.634658\n",
" 0.636347 -0.0144316 -0.771268\n",
" 0.0361865 -0.998166 0.0485336"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"U, S, V = svd(A) #singular decomposition"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "c6fdc076-bfc8-46e9-97c2-5427f4c49395",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 BitMatrix:\n",
" 1 0 0\n",
" 0 1 0\n",
" 0 0 1"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"3×3 BitMatrix:\n",
" 1 0 0\n",
" 0 1 0\n",
" 0 0 1"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"3×3 BitMatrix:\n",
" 0 0 0\n",
" 0 0 0\n",
" 0 0 0"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(U'*U .> 0.00001)\n",
"display(V'*V .> 0.00001)\n",
"display(A - (U * Diagonal(S) * V') .> 0.00001)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "37b7a415-d7ea-4d00-b704-7d79e562d899",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 0.019506810588402184\n",
" 0.2433655332706062\n",
" 0.9153216003851385"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 1.1253009779161582\n",
" 0.27878010490969024\n",
" 0.013851237498536341"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(eigvals(A))\n",
"display(svd(A).S)\n",
"# singular values are more spread out"
]
},
{
"cell_type": "code",
"execution_count": 39,
"id": "5a4893c6-074c-4479-ac8a-9ad8db295d0a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"((3, 3), (3,), (5, 3))"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"((5, 3), (3,), (3, 3))"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# svd exists also for rectangular matrices\n",
"U, S, V = svd(rand(3, 5))\n",
"display((size(U), size(S), size(V)))\n",
"U, S, V = svd(rand(5, 3))\n",
"display((size(U), size(S), size(V)))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "380565b8-f15d-4f3f-9a0f-85fc82809bdd",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 2.302775637731995\n",
" 1.3027756377319948\n",
" 0.0"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"svdvals([2 0 1;0 1 1; 0 0 0; 0 0 0]) # rank is 2 -> so 2 non zero values"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "deb05175-6a64-4c40-ab5c-e614ef85797f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}\n",
"values:\n",
"3-element Vector{Float64}:\n",
" 4.440892098500626e-15\n",
" 1.697224362268007\n",
" 5.302775637731994\n",
"vectors:\n",
"3×3 Matrix{Float64}:\n",
" -0.333333 0.444872 -0.831251\n",
" -0.666667 -0.734656 -0.125841\n",
" 0.666667 -0.51222 -0.541467"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"3×3 Matrix{Float64}:\n",
" 4.0 0.0 2.0\n",
" 0.0 1.0 1.0\n",
" 2.0 1.0 2.0"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"A = [2 0 1;0 1 1; 0 0 0];\n",
"display(eigen(A' * A))\n",
"U, S, V = svd(A);\n",
"display(V * Diagonal(S)' * Diagonal(S) * V')"
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "60516aff-16aa-41e6-a56d-6cad04547e55",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×3 Matrix{Float64}:\n",
" 1.0 -1.0 0.0\n",
" 0.0 1.0 0.0\n",
" 0.0 0.0 0.5"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"3×3 Matrix{Float64}:\n",
" 1.0 -1.0 0.0\n",
" -1.11022e-16 1.0 0.0\n",
" 0.0 0.0 0.5"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"A = [1 1 0; 0 1 0; 0 0 2];\n",
"display(inv(A))\n",
"F = svd(A);\n",
"display(F.V * Diagonal(F.S .^ -1) * F.U') # we can get the inverse from the svd"
]
},
{
"cell_type": "code",
"execution_count": 46,
"id": "cd9e606f-957d-4803-9ab8-f7b23c455e13",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(2.0, false)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"\"frobenius norm: 2.6457513110645907\""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display((opnorm(A), opnorm(A) == norm(A))) # use opnorm for the matrix norm\n",
"display(\"frobenius norm: $(norm(A))\")"
]
},
{
"cell_type": "code",
"execution_count": 57,
"id": "9c469ac2-4d2e-4eba-9c6e-c63ec1fc06ae",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Matrix{Float64}:\n",
" 1.27357 1.80721\n",
" 2.87898 4.08529"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"1"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Eckhart-Young theorem\n",
"# the min X with respect to opnorm(A-X) with rank less than or equal to k is:\n",
"A = [1 2;3 4];\n",
"F = svd(A);\n",
"X = F.U[:,1] .* F.S[1, 1] * F.V[:, 1]';\n",
"display(X)\n",
"display(rank(X))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.9.3",
"language": "julia",
"name": "julia-1.9"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.9.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}

67
10-04/optQf.m Normal file
View File

@ -0,0 +1,67 @@
function v = optQf( Q , q )
%function v = optQf( Q , q )
%
% Given the data structure encoding a Quadratic function
%
% f( x ) = (1/2) x^T * Q * x + q * x
%
% returns the value of f() in the "putative minimum" (or maximum), i.e.,
% the point
%
% xStar = Q \ -q;
%
% Input:
%
% - Q ([ n x n ] real symmetric matrix, not necessarily positive
% semidefinite): the quadratic part of f
%
% - q ([ n x 1 ] real column vector): the linear part of f
%
% Output:
%
% - f( xStar )
%
%{
=======================================
Author: Antonio Frangioni
Date: 30-09-22
Version 0.10
Copyright Antonio Frangioni
=======================================
%}
if ~ isreal( Q )
error( 'Q not a real matrix' );
end
n = size( Q , 1 );
if n <= 1
error( 'Q is too small' );
end
if n ~= size( Q , 2 )
error( 'Q is not square' );
end
if ~ isreal( q )
error( 'q not a real vector' );
end
if size( q , 2 ) ~= 1
error( 'q is not a (column) vector' );
end
if size( q , 1 ) ~= n
error( 'q size does not match with Q' );
end
% "solve the problem" - - - - - - - - - - - - - - - - - - - - - - - - - - -
% ignore if Q is singular
xStar = Q \ -q;
v = 0.5 * xStar' * Q * xStar + q' * xStar;
end

11
10-04/plotQ.jl Normal file
View File

@ -0,0 +1,11 @@
using LinearAlgebra
using Plots
function plotQ(Q::Matrix, q::Vector, xyrange::Tuple{Tuple{Number, Number}, Tuple{Number, Number}})::Plots.Plot
f(x, y) = dot((1/2) * hcat(x, y) * Q, vcat(x, y)) + dot(q', vcat(x, y))
xrange = LinRange(xyrange[1][1], xyrange[2][1], 200)
yrange = LinRange(xyrange[1][2], xyrange[2][2], 200)
z = @. f(xrange', yrange)
plt = contour(xrange, yrange, z)
plt
end

11
10-04/plotQ.m Normal file
View File

@ -0,0 +1,11 @@
function plotQ( Q , q , range )
f = @(x,y) 0.5 * [ x y ] * Q * [ x ; y ] + q' * [ x ; y ];
warning( 'off' , 'all' );
fcontour( f , range , 'LineColor' , 'k' , 'LineWidth' , 1 );
warning( 'on' , 'all' );
end